/*
  cfdtd3d.cpp

  General purpose 3D Electromagnetic Wave Simulation

  Copyright (c) 2011, Bin Guo
  All rights reserved.

  Redistribution and use in source and binary forms,
  with or without modification, are permitted provided
  that the following conditions are met:

    * Redistributions of source code must retain the
	  above copyright notice, this list of conditions
	  and the following disclaimer.
    * Redistributions in binary form must reproduce
	  the above copyright notice, this list of conditions
	  and the following disclaimer in the documentation
	  and/or other materials provided with the distribution.

  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
  CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
  INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
  MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
  STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
  ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


  CHANGES:

  DATE dd/mm/yy          CHANGES
  ----------------------------------------------


 */


#include <stdlib.h>
#include <cstdlib>
#include <stdio.h>
#include <cstdio>
#include <math.h>
#include <exception>
#include <iostream>
#include <time.h>
#include <cstring>

#include "cfdtd3d.h"

using namespace std;


int main(int argc, char *argv[])
{

//All the parameters the code needed during it's running
//is included in the input model file (*.imd)
	if(argc!=3)
	{
		cout<<"Please use the code as:"<<endl;
		cout<<"  cfdtd3d inputfinlename runtimesteps"<<endl;
		return 1;
	}

//Read Data and Initialize
	//Open the input file
    if((fInput=fopen(argv[1],"r"))==NULL)
    {
        cout<<"Open Model File Error!"<<endl;
        return 1;
    }

    //if read data and initialize unsucessfully,
	//show error and exit the running.
    if(ReadData_Init())
    {
		cout<<"Read Data Error!"<<endl;
        return 1;
    }


    T=0;
    StepStart=1;

//input how many time step you want the code calculate.
    nsteps=atoi(argv[2]);

//run the fdtd code
    MainFDTD();

//after running, clear the resource and exit to system.
	DelData();

	return 0;

}




//------------------------------------------
//Read data and initialize subroutine
int ReadData_Init()
{
    int iNumCount=0;
    char StrTemp[201];
    string strReadFileErr="Read Input File Error! Line: ";

//Read Grid Section
//The following 11 lines is a subsection to check is it the right line
//be readed in the input file.
//If not, show error and the number of line.
//Normally, the error happens in the above line.
//You can find this subsection many times at the following code.
    ReadDataFlag=0;
    fgets(StrTemp, 200, fInput);
    iNumCount++;

    //If the code runs in Linux, but the input files are in Windows format,
    //show the error information and return.
    string strTemp1=string(StrTemp);
    string strTemp2=strTemp1.substr(strTemp1.length()-2,strTemp1.length()-1);
    if(!strcmp(strTemp2.c_str(),"\r\n"))
    {
        cout<<"The .imd or .grid file is in Windows format,"<<endl;
        cout<<"please convert it to Linux format, and run the code again."<<endl;
        return 1;
    }

    if(strcmp(StrTemp,"% grid section\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"MaxX, MaxY, MaxZ, dx, dy, dz\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }

//read the maxima X Y Z, and dx, dy, dz
//and find the minimu volue around dx/dy/dz.
    fscanf(fInput,"%d, %d, %d, %lf, %lf, %lf\n", &MaxX,&MaxY,&MaxZ, &dx,&dy,&dz);
    iNumCount++;
    min_dxyz=dx;
	if(min_dxyz>dy) min_dxyz=dy;
	if(min_dxyz>dz) min_dxyz=dz;
    dt=min_dxyz/(2.*C0);
    dx_rate=min_dxyz/dx;
    dy_rate=min_dxyz/dy;
    dz_rate=min_dxyz/dz;

//Read Sourece Section
//Check lines
    iNumCount++;
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"% source section\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"(SourceX1,SourceY1,SourceZ1) (SourceX2,SourceY2,SourceZ2) SourceVector(1x/2y/3z)\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"SourceType, Hard(0)/Soft(1)/Plane(2)/Transfer(3)\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"Source Parameters\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }

//read the position of source
    fscanf(fInput,"(%d,%d,%d) (%d,%d,%d) %d\n",&SourceX1,&SourceY1,&SourceZ1,
        &SourceX2,&SourceY2,&SourceZ2,&SourceVector);
    iNumCount++;

//check the position parameters and find the normal direction of the source port.
    if((SourceX1!=SourceX2)&&(SourceY1!=SourceY2)&&(SourceZ1==SourceZ2))
    {
        SourceNormal=3;
        if(SourceNormal==SourceVector)
        {
            cout<<"SourceNormal==SourceVector"<<endl;
            return 1;
        }
    }
    else if((SourceX1==SourceX2)&&(SourceY1!=SourceY2)&&(SourceZ1!=SourceZ2))
    {
        SourceNormal=1;
        if(SourceNormal==SourceVector)
        {
            cout<<"SourceNormal==SourceVector"<<endl;
            return 1;
        }
    }
    else if((SourceX1!=SourceX2)&&(SourceY1==SourceY2)&&(SourceZ1!=SourceZ2))
    {
        SourceNormal=2;
        if(SourceNormal==SourceVector)
        {
            cout<<"SourceNormal==SourceVector"<<endl;
            return 1;
        }
    }
    else
    {
        SourceNormal=0;
    }
    fscanf(fInput,"%d, %d\n", &SourceType, &SourceHS);
    if((SourceType==2)&&(SourceNormal==0))
    {
        cout<<"PlaneSource can't be used with point and line source1"<<endl;
        return 1;
    }

//read the source wave form parameters based on the source type.
    iNumCount++;
    switch(SourceType)
    {
        case 1:
            fscanf(fInput,"%lf, %lf, %lf\n", &UWB1tao,&UWB1delta,&UWB1t0);
            iNumCount++;
        	UWB1tr=UWB1tao/4.5;
        	UWB1tp=UWB1delta/0.69-UWB1tr;
            break;

        case 2:
            fscanf(fInput,"%lf, %lf\n", &UWB2alpha,&UWB2beita);
    	    iNumCount++;
            UWB2t0=log(UWB2beita/UWB2alpha)/(UWB2beita-UWB2alpha);
	        UWB2E0=UWB2K/(exp(-UWB2alpha*UWB2t0)-exp(-UWB2beita*UWB2t0));
            break;

        case 3:
            fscanf(fInput,"%lf, %lf\n", &GaussSpread, &Gausst0);
            iNumCount++;
            GaussSpread=GaussSpread*1e-12;
            Gausst0=Gausst0*1e-12;
            break;

        case 4:
            fscanf(fInput, "%lf\n", &SinFreq);
            iNumCount++;
            SinFreq=SinFreq*1.E9;
            break;

        case 5:
            fscanf(fInput, "%lf, %lf, %lf\n", &GaussSpread,&Gausst0,&SinFreq);
            iNumCount++;
            GaussSpread=GaussSpread*1e-12;
            Gausst0=Gausst0*1e-12;
            SinFreq=SinFreq*1.E9;
            break;

         case 6:
            fscanf(fInput, "%lf\n", &SinFreq);
            iNumCount++;
            SinFreq=SinFreq*1.E6;
            break;

        case 7:
            fscanf(fInput, "%lf, %lf, %lf\n", &RecSint1,&RecSint2,&RecSinf0);
            iNumCount++;
            RecSint1=RecSint1*1e-12;
            RecSint2=RecSint2*1e-12;
            RecSinf0=RecSinf0*1e9;
            break;

       default:
            cout<<"Source Type Error!"<<endl;
            return 1;
    }
    if((SourceHS!=0)&&(SourceHS!=1)&&(SourceHS!=2)&&(SourceHS!=3))
    {
        cout<<"Source Hard/Soft Error!"<<endl;
        return 1;
    }

//Open source.dat, it is a output file, and the source value at
//every time step will be saved in here.
    if((fSource=fopen("Source.dat","w"))==NULL)
    {
        cout<<"Open Sorce.dat error!"<<endl;
        return 1;
    }

//if the source is a plane type sourec
//general the assistant array.
    if(SourceHS==2)
        switch(SourceNormal)
        {
            case 1:
                E_inc=new double[MaxX];
                H_inc=new double[MaxX];

                for(i=0;i<MaxX;i++)
                {
                    E_inc[i]=0.0;
                    H_inc[i]=0.0;
                }
                break;

            case 2:
                E_inc=new double[MaxY];
                H_inc=new double[MaxY];

                for(j=0;j<MaxY;j++)
                {
                    E_inc[j]=0.0;
                    H_inc[j]=0.0;
                }
                break;

            case 3:
                E_inc=new double[MaxZ];
                H_inc=new double[MaxZ];

                for(k=0;k<MaxZ;k++)
                {
                    E_inc[k]=0.0;
                    H_inc[k]=0.0;
                }
                break;
        }

//Read Output Point Section
//Check lines
    iNumCount++;
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"% output section\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"Num_Output_Point(>=2)\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"(Point1X,Point1Y,Point1Z)\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"(Point2X,Point2Y,Point2Z)\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"...\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fscanf(fInput,"%d\n", &NumWatch);
    iNumCount++;

//If the input file show that the user want to save some
//field information in a point, then open a file name
//watch1.dat and save the value at time step.
//only two points can be saved.
    if(NumWatch>0)
    {
        WatchX=new int[NumWatch];
        WatchY=new int[NumWatch];
        WatchZ=new int[NumWatch];

        for(i=0;i<NumWatch;i++)
		{
            fscanf(fInput,"(%d,%d,%d)\n", &WatchX[i], &WatchY[i], &WatchZ[i]);
			iNumCount++;
		}
    }

    if((fWatch=fopen("Output_Points.dat","w"))==NULL)
    {
        cout<<"Open Output_Point file error!"<<endl;
        return 1;
    }

//Read Watching Port Section
//check lines
    iNumCount++;
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"Num_Output_Port(>=0) both_Voltage_and_Current\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"(PortVX1,PortVY1,PortVZ1)(PortVX2,PortVY2,PortVZ2) (PortIX1,PortIY1,PortIZ1)(PortIX2,PortIY2,PortIZ2)\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"(PortVX1,PortVY1,PortVZ1)(PortVX2,PortVY2,PortVZ2) (PortIX1,PortIY1,PortIZ1)(PortIX2,PortIY2,PortIZ2)\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"......\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }

//read the port information where the user want to save the
//electric voltage and electric current.
//the values will be saved in Vport.dat and Iport1.dat Iport2.dat
//In these files, the 1st coulmn is the time steps and the other
//coulmns are the values at every port.
//The electric current in Iport1.dat is calculated using the magnetic
//field after half-space-step of the electric field which are used to
//calculated the electric voltage.
//The electric current in Iport2.dat is calculated using the magnetic
//field before half-space-step of the electric field which are used to
//calculated the electric voltage.
    fscanf(fInput, "%d\n", &NumPort);
    iNumCount++;

    if(NumPort>0)
    {
        PortVX1=new int[NumPort];
        PortVY1=new int[NumPort];
        PortVZ1=new int[NumPort];
        PortVX2=new int[NumPort];
        PortVY2=new int[NumPort];
        PortVZ2=new int[NumPort];

        PortIX1=new int[NumPort];
        PortIY1=new int[NumPort];
        PortIZ1=new int[NumPort];
        PortIX2=new int[NumPort];
        PortIY2=new int[NumPort];
        PortIZ2=new int[NumPort];


        for(i=0;i<NumPort;i++)
        {
            fscanf(fInput, "(%d,%d,%d)(%d,%d,%d) (%d,%d,%d)(%d,%d,%d)\n",
                &PortVX1[i], &PortVY1[i], &PortVZ1[i],
                &PortVX2[i], &PortVY2[i], &PortVZ2[i],
                &PortIX1[i], &PortIY1[i], &PortIZ1[i],
                &PortIX2[i], &PortIY2[i], &PortIZ2[i]);
            iNumCount++;
        }

        if((fVPort=fopen("VPort.dat","w"))==NULL)
        {
            cout<<"Open VPort File error!"<<endl;
            return 1;
        }
        if((fIPort1=fopen("IPort1.dat","w"))==NULL)
        {
            cout<<"Open IPort1 File error!"<<endl;
            return 1;
        }
        if((fIPort2=fopen("IPort2.dat","w"))==NULL)
        {
            cout<<"Open IPort2 File error!"<<endl;
            return 1;
        }
    }

//Read Watching Plane Section
//check lines
    iNumCount++;
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"Num_Output_Field_on_a_Plane\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"Normal(1x:yz plane / 2y:zx plane /3z:xy plane), Hight, TimeStep\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"......\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }

//read the plane information where the user want to save the
//electric field in this plane at the appointed time step.
//The electric field in every plane will be saved in a file
//named as fieldn.dat. where n is a number as 1,2,3,...
//The number of plane will not be greater than 20 because of
//the limitation of Windows
    char chFieldPlaneIndex[3];
    string strFieldPlaneIndex;
    string strFieldPlaneName;

    fscanf(fInput, "%d\n", &NumFieldPlane);
    iNumCount++;

    if(NumFieldPlane>0)
    {
        FieldPlaneNormal=new int [NumFieldPlane];
        FieldPlaneVector=new int [NumFieldPlane];
        FieldPlaneTimeStep=new int [NumFieldPlane];
        fFieldPlane=new FILE *[NumFieldPlane];

        for(i=0;i<NumFieldPlane;i++)
        {
            fscanf(fInput,"%d, %d, %d\n",&FieldPlaneNormal[i],
                &FieldPlaneVector[i],&FieldPlaneTimeStep[i]);
            iNumCount++;
        }

        for(i=0;i<NumFieldPlane;i++)
        {
            sprintf(chFieldPlaneIndex,"%d",(i+1));
            strFieldPlaneIndex=string(chFieldPlaneIndex);
            strFieldPlaneName="Field"+strFieldPlaneIndex+".dat";

            if((fFieldPlane[i]=fopen(strFieldPlaneName.c_str(),"w"))==NULL)
            {
                cout<<"Open FieldPlane File error!"<<endl;
                return 1;
            }
        }
    }

//Read Medium Parameter Section
//Initialize the medium paramenter array
//The medium type will not be greater than 500
//in a case
    for(i=0;i<MaxMedium;i++)
    {
        Medium_EpsilonR[i]=1.0;
		Medium_EpsilonS[i]=1.0;
        Medium_MiuR[i]=1.0;
        Medium_Sigma[i]=0.0;
		Medium_DebyeTao[i]=-1.0;
        Medium_Conform[i]=0;
        Medium_ConType[i]=1;
    }

//check lines
    iNumCount++;
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"% media section\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"NumMedium\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"MediumID, Epsilon_r, Epsilon_s, Miu_r, Sigma(-1:Infinite), DebyeTao\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"......\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }

//read the information of the medium
    fscanf(fInput,"%d\n", &NumMedium);
    iNumCount++;

    if(NumMedium>0)
    {
        for(i=0;i<NumMedium;i++)
        {
            fscanf(fInput,"%d",&MediumID);
            fscanf(fInput,", %lf, %lf, %lf, %lf, %lf, %d, %d\n",
                &Medium_EpsilonR[MediumID], &Medium_EpsilonS[MediumID],
                &Medium_MiuR[MediumID], &Medium_Sigma[MediumID],
                &Medium_DebyeTao[MediumID],
                &Medium_Conform[MediumID], &Medium_ConType[MediumID]);
            Medium_DebyeTao[MediumID]=Medium_DebyeTao[MediumID]*1E-12;
            iNumCount++;
        }
    }

//Read PEC Section
//Read Normal to X Plane Information
    iNumCount++;
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"% PEC section\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"NumNormalX (PEC on yz plane)\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"xi, y1, y2, z1, z2\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"......\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fscanf(fInput,"%d\n", &NumNormalX);
    iNumCount++;

    if(NumNormalX>0)
    {
        try
        {
            NormalX=new int *[NumNormalX];
            for(i=0;i<NumNormalX;i++) NormalX[i]=new int[5];
        }
        catch (std::bad_alloc)
        {
            cout<<"NormalX could not be allocated."<<endl;
            return(0);
        }

        for(i=0;i<NumNormalX;i++)
        {
            fscanf(fInput,"%d, %d, %d, %d, %d\n",
                &NormalX[i][0],&NormalX[i][1],&NormalX[i][2],
                &NormalX[i][3],&NormalX[i][4]);
            iNumCount++;
        }
    }

//Read Normal to Y Plane Information
    iNumCount++;
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"NumNormalY (PEC on xz plane)\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"yi, x1, x2, z1, z2\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"......\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fscanf(fInput,"%d\n", &NumNormalY);
    iNumCount++;

    if(NumNormalY>0)
    {
        try
        {
            NormalY=new int *[NumNormalY];
            for(i=0;i<NumNormalY;i++) NormalY[i]=new int[5];
        }
        catch (std::bad_alloc)
        {
            cout<<"NormalY could not be allocated."<<endl;
            return 1;
        }

        for(i=0;i<NumNormalY;i++)
        {
            fscanf(fInput,"%d, %d, %d, %d, %d\n",
                &NormalY[i][0],&NormalY[i][1],&NormalY[i][2],
                &NormalY[i][3],&NormalY[i][4]);
            iNumCount++;
        }
    }

//Read Normal to Z Plane Information
    iNumCount++;
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"NumNormalZ (PEC on xy plane)\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"zi, x1, x2, y1, y2\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"......\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fscanf(fInput," %d\n", &NumNormalZ);
    iNumCount++;

    if(NumNormalZ>0)
    {
        try
        {
            NormalZ=new int *[NumNormalZ];
            for(i=0;i<NumNormalZ;i++) NormalZ[i]=new int[5];
        }
        catch (std::bad_alloc)
        {
            cout<<"NormalZ could not be allocated."<<endl;
            return(0);
        }

        for(i=0;i<NumNormalZ;i++)
        {
            fscanf(fInput,"%d, %d, %d, %d, %d\n",
                &NormalZ[i][0],&NormalZ[i][1],&NormalZ[i][2],
                &NormalZ[i][3],&NormalZ[i][4]);
            iNumCount++;
        }
    }

//Read PML Section
//Initialize the arrays will be used to updata PML
    iNumCount++;
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"% PML section\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"PML_lx, PML_ly, PML_lz\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fscanf(fInput,"%d, %d, %d\n",&PML_lx, &PML_ly, &PML_lz);
    iNumCount++;

//Read Model Section
    iNumCount++;
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"% Model section\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    fgets(StrTemp, 200, fInput);
    iNumCount++;
    if(strcmp(StrTemp,"FileName\n"))
    {
        cout<<strReadFileErr<<iNumCount<<endl;
        return 1;
    }
    char chModelName[200];
    fgets(chModelName, 200, fInput);
    string strNameTemp=string(chModelName);
    string strModelName=strNameTemp.substr(0,strNameTemp.length()-1);

//Open Model file to read the Meshed Model
    if((fModel=fopen(strModelName.c_str(),"r"))==NULL)
    {
        cout<<"Open Model file error!"<<endl;
        return 1;
    }
    int mod_maxx, mod_maxy, mod_maxz;
    fscanf(fModel,"%d, %d, %d\n", &mod_maxx, &mod_maxy,&mod_maxz);
    if((mod_maxx!=MaxX)|(mod_maxy!=MaxY)|(mod_maxz!=MaxZ))
    {
        cout<<"The size of the Model file does not match the Input file!"<<endl;
        return 1;
    }


//Init Model array based on the maximum X, Y, Z
    try
    {
        Model=new int **[MaxX];
        for(i=0;i<MaxX;i++) Model[i]=new int *[MaxY];
        for(i=0;i<MaxX;i++)
            for(j=0;j<MaxY;j++)
                Model[i][j]=new int[MaxZ];
    }
    catch (std::bad_alloc)
    {
        cout<<"Model could not be allocated."<<endl;
        return 1;
    }

//Read Model Data to the model array
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            for(k=0;k<MaxZ;k++)
                Model[i][j][k]=0;

    for(k=0;k<MaxZ;k++)
        for(j=0;j<MaxY;j++)
        {
            for(i=0;i<MaxX;i++)
                fscanf(fModel,"%d",&Model[i][j][k]);
            fscanf(fModel,"\n");
        }

    fclose(fInput);
    fclose(fModel);
//Read Data End


//Init Dx, Dy, Dz array
        try
        {
            Dx=new double **[MaxX];
            for(i=0;i<MaxX;i++) Dx[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Dx[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Dx could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Dy=new double **[MaxX];
            for(i=0;i<MaxX;i++) Dy[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Dy[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Dy could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Dz=new double **[MaxX];
            for(i=0;i<MaxX;i++) Dz[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Dz[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Dz could not be allocated."<<endl;
            return 1;
        }

//Init Ex, Ey, Ez array
        try
        {
            Ex=new double **[MaxX];
            for(i=0;i<MaxX;i++) Ex[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Ex[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Ex could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Ey=new double **[MaxX];
            for(i=0;i<MaxX;i++) Ey[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Ey[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Ey could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Ez=new double **[MaxX];
            for(i=0;i<MaxX;i++) Ez[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Ez[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Ez could not be allocated."<<endl;
            return 1;
        }

//Init gax, gay, gaz array
        try
        {
            gax=new double **[MaxX];
            for(i=0;i<MaxX;i++) gax[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    gax[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"gax could not be allocated."<<endl;
            return 1;
        }
        try
        {
            gay=new double **[MaxX];
            for(i=0;i<MaxX;i++) gay[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    gay[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"gay could not be allocated."<<endl;
            return 1;
        }
        try
        {
            gaz=new double **[MaxX];
            for(i=0;i<MaxX;i++) gaz[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    gaz[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"gaz could not be allocated."<<endl;
            return 1;
        }

//Init gbx, gby, gbz array
        try
        {
            gbx=new double **[MaxX];
            for(i=0;i<MaxX;i++) gbx[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    gbx[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"gbx could not be allocated."<<endl;
            return 1;
        }
        try
        {
            gby=new double **[MaxX];
            for(i=0;i<MaxX;i++) gby[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    gby[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"gby could not be allocated."<<endl;
            return 1;
        }
        try
        {
            gbz=new double **[MaxX];
            for(i=0;i<MaxX;i++) gbz[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    gbz[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"gbz could not be allocated."<<endl;
            return 1;
        }

//Init gcx, gcy, gcz array
        try
        {
            gcx=new double **[MaxX];
            for(i=0;i<MaxX;i++) gcx[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    gcx[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"gcx could not be allocated."<<endl;
            return 1;
        }
        try
        {
            gcy=new double **[MaxX];
            for(i=0;i<MaxX;i++) gcy[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    gcy[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"gcy could not be allocated."<<endl;
            return 1;
        }
        try
        {
            gcz=new double **[MaxX];
            for(i=0;i<MaxX;i++) gcz[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    gcz[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"gcz could not be allocated."<<endl;
            return 1;
        }

//Init gbx, gby, gbz array
        try
        {
            gdx=new double **[MaxX];
            for(i=0;i<MaxX;i++) gdx[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    gdx[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"gdx could not be allocated."<<endl;
            return 1;
        }
        try
        {
            gdy=new double **[MaxX];
            for(i=0;i<MaxX;i++) gdy[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    gdy[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"gdy could not be allocated."<<endl;
            return 1;
        }
        try
        {
            gdz=new double **[MaxX];
            for(i=0;i<MaxX;i++) gdz[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    gdz[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"gdz could not be allocated."<<endl;
            return 1;
        }

//Init Ix, Iy, Iz array
        try
        {
            Ix=new double **[MaxX];
            for(i=0;i<MaxX;i++) Ix[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Ix[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Ix could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Iy=new double **[MaxX];
            for(i=0;i<MaxX;i++) Iy[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Iy[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Iy could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Iz=new double **[MaxX];
            for(i=0;i<MaxX;i++) Iz[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Iz[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Iz could not be allocated."<<endl;
            return 1;
        }

//Init Sx, Sy, Sz array
        try
        {
            Sx=new double **[MaxX];
            for(i=0;i<MaxX;i++) Sx[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Sx[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Sx could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Sy=new double **[MaxX];
            for(i=0;i<MaxX;i++) Sy[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Sy[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Sy could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Sz=new double **[MaxX];
            for(i=0;i<MaxX;i++) Sz[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Sz[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Sz could not be allocated."<<endl;
            return 1;
        }

//Init Hx, Hy, Hz array
        try
        {
            Hx=new double **[MaxX];
            for(i=0;i<MaxX;i++) Hx[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Hx[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Hx could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Hy=new double **[MaxX];
            for(i=0;i<MaxX;i++) Hy[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Hy[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Hy could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Hz=new double **[MaxX];
            for(i=0;i<MaxX;i++) Hz[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Hz[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Hz could not be allocated."<<endl;
            return 1;
        }

//Init HB array
        try
        {
            HB=new double **[MaxX];
            for(i=0;i<MaxX;i++) HB[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    HB[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"HB could not be allocated."<<endl;
            return 1;
        }


//Init Id array
        try
        {
            Idxl=new double **[PML_lx];
            for(i=0;i<PML_lx;i++) Idxl[i]=new double *[MaxY];
            for(i=0;i<PML_lx;i++)
                for(j=0;j<MaxY;j++)
                    Idxl[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Idxl could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Idxh=new double **[PML_lx];
            for(i=0;i<PML_lx;i++) Idxh[i]=new double *[MaxY];
            for(i=0;i<PML_lx;i++)
                for(j=0;j<MaxY;j++)
                    Idxh[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Idxh could not be allocated."<<endl;
            return 1;
        }

        try
        {
            Idyl=new double **[MaxX];
            for(i=0;i<MaxX;i++) Idyl[i]=new double *[PML_ly];
            for(i=0;i<MaxX;i++)
                for(j=0;j<PML_ly;j++)
                    Idyl[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Idyl could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Idyh=new double **[MaxX];
            for(i=0;i<MaxX;i++) Idyh[i]=new double *[PML_ly];
            for(i=0;i<MaxX;i++)
                for(j=0;j<PML_ly;j++)
                    Idyh[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Idyh could not be allocated."<<endl;
            return 1;
        }

        try
        {
            Idzl=new double **[MaxX];
            for(i=0;i<MaxX;i++) Idzl[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Idzl[i][j]=new double[PML_lz];
        }
        catch (std::bad_alloc)
        {
            cout<<"Idzl could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Idzh=new double **[MaxX];
            for(i=0;i<MaxX;i++) Idzh[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Idzh[i][j]=new double[PML_lz];
        }
        catch (std::bad_alloc)
        {
            cout<<"Idzh could not be allocated."<<endl;
            return 1;
        }

//Init Ih array
        try
        {
            Ihxl=new double **[PML_lx];
            for(i=0;i<PML_lx;i++) Ihxl[i]=new double *[MaxY];
            for(i=0;i<PML_lx;i++)
                for(j=0;j<MaxY;j++)
                    Ihxl[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Ihxl could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Ihxh=new double **[PML_lx];
            for(i=0;i<PML_lx;i++) Ihxh[i]=new double *[MaxY];
            for(i=0;i<PML_lx;i++)
                for(j=0;j<MaxY;j++)
                    Ihxh[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Ihxh could not be allocated."<<endl;
            return 1;
        }

        try
        {
            Ihyl=new double **[MaxX];
            for(i=0;i<MaxX;i++) Ihyl[i]=new double *[PML_ly];
            for(i=0;i<MaxX;i++)
                for(j=0;j<PML_ly;j++)
                    Ihyl[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Ihyl could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Ihyh=new double **[MaxX];
            for(i=0;i<MaxX;i++) Ihyh[i]=new double *[PML_ly];
            for(i=0;i<MaxX;i++)
                for(j=0;j<PML_ly;j++)
                    Ihyh[i][j]=new double[MaxZ];
        }
        catch (std::bad_alloc)
        {
            cout<<"Ihyh could not be allocated."<<endl;
            return 1;
        }

        try
        {
            Ihzl=new double **[MaxX];
            for(i=0;i<MaxX;i++) Ihzl[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Ihzl[i][j]=new double[PML_lz];
        }
        catch (std::bad_alloc)
        {
            cout<<"Ihzl could not be allocated."<<endl;
            return 1;
        }
        try
        {
            Ihzh=new double **[MaxX];
            for(i=0;i<MaxX;i++) Ihzh[i]=new double *[MaxY];
            for(i=0;i<MaxX;i++)
                for(j=0;j<MaxY;j++)
                    Ihzh[i][j]=new double[PML_lz];
        }
        catch (std::bad_alloc)
        {
            cout<<"Ihzh could not be allocated."<<endl;
            return 1;
        }


//Initialize Matrices
for(i=0;i<MaxX;i++)
    for(j=0;j<MaxY;j++)
		for(k=0;k<MaxZ;k++)
        {
            Dx[i][j][k]=0.0;
            Dy[i][j][k]=0.0;
            Dz[i][j][k]=0.0;
            Ex[i][j][k]=0.0;
            Ey[i][j][k]=0.0;
            Ez[i][j][k]=0.0;
            gax[i][j][k]=1.0;
            gay[i][j][k]=1.0;
            gaz[i][j][k]=1.0;
            gbx[i][j][k]=0.0;
            gby[i][j][k]=0.0;
            gbz[i][j][k]=0.0;
            gcx[i][j][k]=0.0;
            gcy[i][j][k]=0.0;
            gcz[i][j][k]=0.0;
            gdx[i][j][k]=0.0;
            gdy[i][j][k]=0.0;
            gdz[i][j][k]=0.0;
            Ix[i][j][k]=0.0;
            Iy[i][j][k]=0.0;
            Iz[i][j][k]=0.0;
            Sx[i][j][k]=0.0;
            Sy[i][j][k]=0.0;
            Sz[i][j][k]=0.0;
            Hx[i][j][k]=0.0;
            Hy[i][j][k]=0.0;
            Hz[i][j][k]=0.0;
            HB[i][j][k]=1.0;
        }

for(i=0;i<PML_lx;i++)
    for(j=0;j<MaxY;j++)
        for(k=0;k<MaxZ;k++)
        {
            Idxl[i][j][k]=0.0;
            Idxh[i][j][k]=0.0;
            Ihxl[i][j][k]=0.0;
            Ihxh[i][j][k]=0.0;
        }
for(i=0;i<MaxX;i++)
    for(j=0;j<PML_ly;j++)
        for(k=0;k<MaxZ;k++)
        {
            Idyl[i][j][k]=0.0;
            Idyh[i][j][k]=0.0;
            Ihyl[i][j][k]=0.0;
            Ihyh[i][j][k]=0.0;
        }
for(i=0;i<MaxX;i++)
    for(j=0;j<MaxY;j++)
        for(k=0;k<PML_lz;k++)
        {
            Idzl[i][j][k]=0.0;
            Idzh[i][j][k]=0.0;
            Ihzl[i][j][k]=0.0;
            Ihzh[i][j][k]=0.0;
        }
//Initialize Matrices End

//Generate HB, gax,gay,gaz, gbx,gby,gbz
for(k=0;k<MaxZ;k++)
    for(j=0;j<MaxY;j++)
        for(i=0;i<MaxX;i++)
        {
            Epsilon_r=Medium_EpsilonR[Model[i][j][k]];
            Epsilon_s=Medium_EpsilonS[Model[i][j][k]];
            Miu_r=Medium_MiuR[Model[i][j][k]];
            Sigma=Medium_Sigma[Model[i][j][k]];
			DebyeTao=Medium_DebyeTao[Model[i][j][k]];
            Conform=Medium_Conform[Model[i][j][k]];
            ConType=Medium_ConType[Model[i][j][k]];

            HB[i][j][k]=1./Miu_r;

            if((Conform==0)&&(Sigma>=0.0))
            {
                if(DebyeTao>0.0)
				{
					gdx[i][j][k]=(Epsilon_s-Epsilon_r)*dt/DebyeTao;
					gdy[i][j][k]=(Epsilon_s-Epsilon_r)*dt/DebyeTao;
					gdz[i][j][k]=(Epsilon_s-Epsilon_r)*dt/DebyeTao;
				}
				else
				{
					gdx[i][j][k]=0.0;
					gdy[i][j][k]=0.0;
					gdz[i][j][k]=0.0;
				}

                gax[i][j][k]=1.0/(Epsilon_r+(Sigma*dt/Eps0)+gdx[i][j][k]);
                gay[i][j][k]=1.0/(Epsilon_r+(Sigma*dt/Eps0)+gdy[i][j][k]);
                gaz[i][j][k]=1.0/(Epsilon_r+(Sigma*dt/Eps0)+gdz[i][j][k]);

                gbx[i][j][k]=Sigma*dt/Eps0;
                gby[i][j][k]=Sigma*dt/Eps0;
                gbz[i][j][k]=Sigma*dt/Eps0;

                if(DebyeTao>0.0)
				{
					gcx[i][j][k]=exp(-1.*dt/DebyeTao);
					gcy[i][j][k]=exp(-1.*dt/DebyeTao);
					gcz[i][j][k]=exp(-1.*dt/DebyeTao);
				}
				else
				{
					gcx[i][j][k]=0.0;
					gcy[i][j][k]=0.0;
					gcz[i][j][k]=0.0;
				}
            }
            else if((Conform==0)&&(Sigma<0.0))
            {
                gax[i][j][k]=0.0;
                gay[i][j][k]=0.0;
                gaz[i][j][k]=0.0;
                gbx[i][j][k]=0.0;
                gby[i][j][k]=0.0;
                gbz[i][j][k]=0.0;
                gcx[i][j][k]=0.0;
                gcy[i][j][k]=0.0;
                gcz[i][j][k]=0.0;
                gdx[i][j][k]=0.0;
                gdy[i][j][k]=0.0;
                gdz[i][j][k]=0.0;
            }
            else if(Conform!=0)
            {
                switch(Conform)
                {
                    case 1:
                        if(ConType==11)
                        {
                            gax[i][j][k]=0.0;
                            gay[i][j][k]=1.0/(Epsilon_r+(Sigma*dt/Eps0));
                            gaz[i][j][k]=0.0;
                            gbx[i][j][k]=0.0;
                            gby[i][j][k]=Sigma*dt/Eps0;
                            gbz[i][j][k]=0.0;
                        }
                        else if(ConType==12)
                        {
                            gax[i][j][k]=1.0/(Epsilon_r+(Sigma*dt/Eps0));
                            gay[i][j][k]=0.0;
                            gaz[i][j][k]=0.0;
                            gbx[i][j][k]=Sigma*dt/Eps0;
                            gby[i][j][k]=0.0;
                            gbz[i][j][k]=0.0;
                        }
                        else if(ConType==21)
                        {
                            gax[i][j][k]=1.0/(Epsilon_r+(Sigma*dt/Eps0));
                            gay[i][j][k]=1.0/(Epsilon_r+(Sigma*dt/Eps0));
                            gaz[i][j][k]=0.0;
                            gbx[i][j][k]=Sigma*dt/Eps0;
                            gby[i][j][k]=Sigma*dt/Eps0;
                            gbz[i][j][k]=0.0;
                        }
                        else if(ConType==22)
                        {
                            gax[i][j][k]=0.0;
                            gay[i][j][k]=0.0;
                            gaz[i][j][k]=0.0;
                            gbx[i][j][k]=0.0;
                            gby[i][j][k]=0.0;
                            gbz[i][j][k]=0.0;
                        }
                    break;

                    case 2:
                        if(ConType==11)
                        {
                            gay[i][j][k]=0.0;
                            gaz[i][j][k]=1.0/(Epsilon_r+(Sigma*dt/Eps0));
                            gax[i][j][k]=0.0;
                            gby[i][j][k]=0.0;
                            gbz[i][j][k]=Sigma*dt/Eps0;
                            gbx[i][j][k]=0.0;
                        }
                        else if(ConType==12)
                        {
                            gay[i][j][k]=1.0/(Epsilon_r+(Sigma*dt/Eps0));
                            gaz[i][j][k]=0.0;
                            gax[i][j][k]=0.0;
                            gby[i][j][k]=Sigma*dt/Eps0;
                            gbz[i][j][k]=0.0;
                            gbx[i][j][k]=0.0;
                        }
                        else if(ConType==21)
                        {
                            gay[i][j][k]=1.0/(Epsilon_r+(Sigma*dt/Eps0));
                            gaz[i][j][k]=1.0/(Epsilon_r+(Sigma*dt/Eps0));
                            gax[i][j][k]=0.0;
                            gby[i][j][k]=Sigma*dt/Eps0;
                            gbz[i][j][k]=Sigma*dt/Eps0;
                            gbx[i][j][k]=0.0;
                        }
                        else if(ConType==22)
                        {
                            gay[i][j][k]=0.0;
                            gaz[i][j][k]=0.0;
                            gax[i][j][k]=0.0;
                            gby[i][j][k]=0.0;
                            gbz[i][j][k]=0.0;
                            gbx[i][j][k]=0.0;
                        }
                    break;

                    case 3:
                        if(ConType==11)
                        {
                            gaz[i][j][k]=0.0;
                            gax[i][j][k]=1.0/(Epsilon_r+(Sigma*dt/Eps0));
                            gay[i][j][k]=0.0;
                            gbz[i][j][k]=0.0;
                            gbx[i][j][k]=Sigma*dt/Eps0;
                            gby[i][j][k]=0.0;
                        }
                        else if(ConType==12)
                        {
                            gaz[i][j][k]=1.0/(Epsilon_r+(Sigma*dt/Eps0));
                            gax[i][j][k]=0.0;
                            gay[i][j][k]=0.0;
                            gbz[i][j][k]=Sigma*dt/Eps0;
                            gbx[i][j][k]=0.0;
                            gby[i][j][k]=0.0;
                        }
                        else if(ConType==21)
                        {
                            gaz[i][j][k]=1.0/(Epsilon_r+(Sigma*dt/Eps0));
                            gax[i][j][k]=1.0/(Epsilon_r+(Sigma*dt/Eps0));
                            gay[i][j][k]=0.0;
                            gbz[i][j][k]=Sigma*dt/Eps0;
                            gbx[i][j][k]=Sigma*dt/Eps0;
                            gby[i][j][k]=0.0;
                        }
                        else if(ConType==22)
                        {
                            gaz[i][j][k]=0.0;
                            gax[i][j][k]=0.0;
                            gay[i][j][k]=0.0;
                            gbz[i][j][k]=0.0;
                            gbx[i][j][k]=0.0;
                            gby[i][j][k]=0.0;
                        }
                    break;

                    default:
                        gax[i][j][k]=0.0;
                        gay[i][j][k]=0.0;
                        gaz[i][j][k]=0.0;
                        gbx[i][j][k]=0.0;
                        gby[i][j][k]=0.0;
                        gbz[i][j][k]=0.0;
                }
            }
        }

//Conductive Plate normal to x/y/z aixl
for(i=0;i<NumNormalX;i++)
{
    for(j=NormalX[i][1];j<=NormalX[i][2];j++)
        for(k=NormalX[i][3];k<=NormalX[i][4]+1;k++)
            gay[NormalX[i][0]][j][k]=0.0;

    for(j=NormalX[i][1];j<=NormalX[i][2]+1;j++)
        for(k=NormalX[i][3];k<=NormalX[i][4];k++)
            gaz[NormalX[i][0]][j][k]=0.0;
}

for(j=0;j<NumNormalY;j++)
{
    for(i=NormalY[j][1];i<=NormalY[j][2];i++)
        for(k=NormalY[j][3];k<=NormalY[j][4]+1;k++)
            gax[i][NormalY[j][0]][k]=0.0;

    for(i=NormalY[j][1];i<=NormalY[j][2]+1;i++)
        for(k=NormalY[j][3];k<=NormalY[j][4];k++)
            gaz[i][NormalY[j][0]][k]=0.0;
}

for(k=0;k<NumNormalZ;k++)
{
    for(i=NormalZ[k][1];i<=NormalZ[k][2];i++)
        for(j=NormalZ[k][3];j<=NormalZ[k][4]+1;j++)
            gax[i][j][NormalZ[k][0]]=0.0;

    for(i=NormalZ[k][1];i<=NormalZ[k][2]+1;i++)
        for(j=NormalZ[k][3];j<=NormalZ[k][4];j++)
            gay[i][j][NormalZ[k][0]]=0.0;
}

// HB, ga, gb, End

//PML Parameter
    //X Direction
    try
    {
        gx3=new double[MaxX];
    }
    catch (std::bad_alloc)
    {
        cout<<"gx3 could not be allocated."<<endl;
        return 1;
    }
    try
    {
        gx2=new double[MaxX];
    }
    catch (std::bad_alloc)
    {
        cout<<"gx2 could not be allocated."<<endl;
        return 1;
    }
    try
    {
        gx1=new double[MaxX];
    }
    catch (std::bad_alloc)
    {
        cout<<"gx1 could not be allocated."<<endl;
        return 1;
    }
    try
    {
        fx1=new double[MaxX];
    }
    catch (std::bad_alloc)
    {
        cout<<"fx1 could not be allocated."<<endl;
        return 1;
    }
    try
    {
        fx2=new double[MaxX];
    }
    catch (std::bad_alloc)
    {
        cout<<"fx2 could not be allocated."<<endl;
        return 1;
    }
    try
    {
        fx3=new double[MaxX];
    }
    catch (std::bad_alloc)
    {
        cout<<"fx3 could not be allocated."<<endl;
        return 1;
    }

    //Y Direction
    try
    {
        gy3=new double[MaxY];
    }
    catch (std::bad_alloc)
    {
        cout<<"gy3 could not be allocated."<<endl;
        return 1;
    }
    try
    {
        gy2=new double[MaxY];
    }
    catch (std::bad_alloc)
    {
        cout<<"gy2 could not be allocated."<<endl;
        return 1;
    }
    try
    {
        gy1=new double[MaxY];
    }
    catch (std::bad_alloc)
    {
        cout<<"gy1 could not be allocated."<<endl;
        return 1;
    }
    try
    {
        fy1=new double[MaxY];
    }
    catch (std::bad_alloc)
    {
        cout<<"fy1 could not be allocated."<<endl;
        return 1;
    }
    try
    {
        fy2=new double[MaxY];
    }
    catch (std::bad_alloc)
    {
        cout<<"fy2 could not be allocated."<<endl;
        return 1;
    }
    try
    {
        fy3=new double[MaxY];
    }
    catch (std::bad_alloc)
    {
        cout<<"fy3 could not be allocated."<<endl;
        return 1;
    }

    //Z Direction
    try
    {
        gz3=new double[MaxZ];
    }
    catch (std::bad_alloc)
    {
        cout<<"gz3 could not be allocated."<<endl;
        return 1;
    }
    try
    {
        gz2=new double[MaxZ];
    }
    catch (std::bad_alloc)
    {
        cout<<"gz2 could not be allocated."<<endl;
        return 1;
    }
    try
    {
        gz1=new double[MaxZ];
    }
    catch (std::bad_alloc)
    {
        cout<<"gz1 could not be allocated."<<endl;
        return 1;
    }
    try
    {
        fz1=new double[MaxZ];
    }
    catch (std::bad_alloc)
    {
        cout<<"fz1 could not be allocated."<<endl;
        return 1;
    }
    try
    {
        fz2=new double[MaxZ];
    }
    catch (std::bad_alloc)
    {
        cout<<"fz2 could not be allocated."<<endl;
        return 1;
    }
    try
    {
        fz3=new double[MaxZ];
    }
    catch (std::bad_alloc)
    {
        cout<<"fz3 could not be allocated."<<endl;
        return 1;
    }

//Init gx,fx value
    for(i=0;i<MaxX;i++)
    {
        gx1[i]=0.0;
        gx2[i]=1.0;
        gx3[i]=1.0;
        fx1[i]=0.0;
        fx2[i]=1.0;
        fx3[i]=1.0;
    }
    for(j=0;j<MaxY;j++)
    {
        gy1[j]=0.0;
        gy2[j]=1.0;
        gy3[j]=1.0;
        fy1[j]=0.0;
        fy2[j]=1.0;
        fy3[j]=1.0;
    }
    for(k=0;k<MaxZ;k++)
    {
        gz1[k]=0.0;
        gz2[k]=1.0;
        gz3[k]=1.0;
        fz1[k]=0.0;
        fz2[k]=1.0;
        fz3[k]=1.0;
    }

    for(i=0;i<PML_lx;i++)
    {
        xnum=PML_lx-i;
        xd=PML_lx;
        xxn=xnum/xd;
        PML_factor=(dt/2./Eps0)*((PML_m+1.)/(150.*Pi*dx*sqrt(PML_epsr)));
        xnpml=PML_factor*pow(xxn,PML_m);
        gx1[i]=xnpml;
        gx1[MaxX-2-i]=xnpml;
        gx2[i]=1.0/(1.0+xnpml);
        gx2[MaxX-1-i]=1.0/(1.0+xnpml);
        gx3[i]=(1.0-xnpml)/(1.0+xnpml);
        gx3[MaxX-1-i]=(1.0-xnpml)/(1.0+xnpml);
        xxn=(xnum-0.5)/xd;
        xnpml=PML_factor*pow(xxn,PML_m);
        fx1[i]=xnpml;
        fx1[MaxX-1-i]=xnpml;
        fx2[i]=1.0/(1.0+xnpml);
        fx2[MaxX-2-i]=1.0/(1.0+xnpml);
        fx3[i]=(1.0-xnpml)/(1.0+xnpml);
        fx3[MaxX-2-i]=(1.0-xnpml)/(1.0+xnpml);
    }

    for(j=0;j<PML_ly;j++)
    {
        ynum=PML_ly-j;
        yd=PML_ly;
        yyn=ynum/yd;
        PML_factor=(dt/2./Eps0)*((PML_m+1.)/(150.*Pi*dy*sqrt(PML_epsr)));
        ynpml=PML_factor*pow(yyn,PML_m);
        gy1[j]=ynpml;
        gy1[MaxY-2-j]=ynpml;
        gy2[j]=1.0/(1.0+ynpml);
        gy2[MaxY-1-j]=1.0/(1.0+ynpml);
        gy3[j]=(1.0-ynpml)/(1.0+ynpml);
        gy3[MaxY-1-j]=(1.0-ynpml)/(1.0+ynpml);
        yyn=(ynum-0.5)/yd;
        ynpml=PML_factor*pow(yyn,PML_m);
        fy1[j]=ynpml;
        fy1[MaxY-1-j]=ynpml;
        fy2[j]=1.0/(1.0+ynpml);
        fy2[MaxY-2-j]=1.0/(1.0+ynpml);
        fy3[j]=(1.0-ynpml)/(1.0+ynpml);
        fy3[MaxY-2-j]=(1.0-ynpml)/(1.0+ynpml);
    }

    for(k=0;k<PML_lz;k++)
    {
        znum=PML_ly-k;
        zd=PML_lz;
        zzn=znum/zd;
        PML_factor=(dt/2./Eps0)*((PML_m+1.)/(150.*Pi*dz*sqrt(PML_epsr)));
        znpml=PML_factor*pow(zzn,PML_m);
        gz1[k]=znpml;
        gz1[MaxZ-2-k]=znpml;
        gz2[k]=1.0/(1.0+znpml);
        gz2[MaxZ-1-k]=1.0/(1.0+znpml);
        gz3[k]=(1.0-znpml)/(1.0+znpml);
        gz3[MaxZ-1-k]=(1.0-znpml)/(1.0+znpml);
        zzn=(znum-0.5)/zd;
        znpml=PML_factor*pow(zzn,PML_m);
        fz1[k]=znpml;
        fz1[MaxZ-1-k]=znpml;
        fz2[k]=1.0/(1.0+znpml);
        fz2[MaxZ-2-k]=1.0/(1.0+znpml);
        fz3[k]=(1.0-znpml)/(1.0+znpml);
        fz3[MaxZ-2-k]=(1.0-znpml)/(1.0+znpml);
    }

    /*FILE *fPMLx,*fPMLy,*fPMLz;
    fPMLx=fopen("PMLx.dat","w");
    fPMLy=fopen("PMLy.dat","w");
    fPMLz=fopen("PMLz.dat","w");
    for(i=0;i<MaxX;i++)
        fprintf(fPMLx,"%16.8E, %16.8E, %16.8E, %16.8E, %16.8E, %16.8E, %16.8E\n",gx1[i]*2.*Eps0/dt,
            gx1[i],gx2[i],gx3[i],fx1[i],fx2[i],fx3[i]);

    for(j=0;j<MaxY;j++)
        fprintf(fPMLy,"%16.8E, %16.8E, %16.8E, %16.8E, %16.8E, %16.8E, %16.8E\n",gy1[j]*2.*Eps0/dt,
            gy1[j],gy2[j],gy3[j],fy1[j],fy2[j],fy3[j]);

    for(k=0;k<MaxZ;k++)
        fprintf(fPMLz,"%16.8E, %16.8E, %16.8E, %16.8E, %16.8E, %16.8E, %16.8E\n",gz1[k]*2.*Eps0/dt,
            gz1[k],gz2[k],gz3[k],fz1[k],fz2[k],fz3[k]);

    fclose(fPMLx);
    fclose(fPMLy);
    fclose(fPMLz);*/

//PML Parameter END

    ReadDataFlag=1;
    OpenFileFlag=0;
    DelDataFlag=0;

    return 0;
}
//------------------------------------------



//------------------------------------------
//Main FDTD
int MainFDTD()
{
    isRunning=1;

//Save the time the code beging to run
    t_start=time(NULL);

//start the main updata
    for(n=StepStart;n<StepStart+nsteps;n++)
    {
		//show the time step
		cout<<"Present step is : "<<n<<" of "<<nsteps<<endl;

        T++;

		//updata D field
        CalculateD();

			//Add source the fdtd domain
            /*if(SourceHS==1)
                AddSourceD();
            else if(SourceHS==2)
                AddSourceE();*/
			if(SourceHS==2)
				AddSourceE();

            BoundD_2_Zero();

            //AddDiodeIdeal();
            //AddDiode1();

		//updata E field
		CalculateE();

            if(SourceHS==0)
                AddSource();
			else if(SourceHS==1)
				AddSoftSource();

            BoundE_2_Zero();

			//let electric field in PEC equrie zero
            ConductorE_2_Zero();

        SavePortData();

		//updata H field
        CalculateH();

            if(SourceHS==2)
                AddSourceH();

            BoundH_2_Zero();


        switch (SourceVector)
        {
            case 1:
             fprintf(fWatch,"%16.8E",dt*T*1e12);
             for(i=0;i<NumWatch;i++)
                fprintf(fWatch," %16.8E",Ex[WatchX[i]][WatchY[i]][WatchZ[i]]);
             fprintf(fWatch,"\n");

             break;

            case 2:
             fprintf(fWatch,"%16.8E",dt*T*1e12);
             for(i=0;i<NumWatch;i++)
                fprintf(fWatch," %16.8E",Ey[WatchX[i]][WatchY[i]][WatchZ[i]]);
             fprintf(fWatch,"\n");

             break;

            case 3:
             fprintf(fWatch,"%16.8E",dt*T*1e12);
             for(i=0;i<NumWatch;i++)
                fprintf(fWatch," %16.8E",Ez[WatchX[i]][WatchY[i]][WatchZ[i]]);
             fprintf(fWatch,"\n");

             break;

            default:
             fprintf(fWatch,"%16.8E",dt*T*1e12);
             for(i=0;i<NumWatch;i++)
                fprintf(fWatch," %16.8E",Ex[WatchX[i]][WatchY[i]][WatchZ[i]]);
             fprintf(fWatch,"\n");
        }


        //Save Field on Plane
        SaveFieldonPlane();

    }

    StepStart=StepStart+nsteps;
    isRunning=0;

//save the end time of the code
//and show how long the code run
	t_end=time(NULL);
	cout<<"Used "<<(t_end-t_start)<<" seconds."<<endl;

    return(0);
}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
//updata D field
void CalculateD()
{
//Calculate D Field
        //Calculate Dx
        for(i=0;i<PML_lx;i++)
            for(j=1;j<MaxY;j++)
                for(k=1;k<MaxZ;k++)
                {
                    curl_h=dy_rate*(Hz[i][j][k]-Hz[i][j-1][k])-
                           dz_rate*(Hy[i][j][k]-Hy[i][j][k-1]);
                    Idxl[i][j][k]=Idxl[i][j][k]+curl_h;
                    Dx[i][j][k]=gy3[j]*gz3[k]*Dx[i][j][k]+gy2[j]*gz2[k]*
                        0.5*(curl_h+gx1[i]*Idxl[i][j][k]);
                }

        for(i=PML_lx;i<MaxX-PML_lx-1;i++)
            for(j=1;j<MaxY;j++)
                for(k=1;k<MaxZ;k++)
                {
                    curl_h=dy_rate*(Hz[i][j][k]-Hz[i][j-1][k])-
                           dz_rate*(Hy[i][j][k]-Hy[i][j][k-1]);
                    Dx[i][j][k]=gy3[j]*gz3[k]*Dx[i][j][k]+
                                gy2[j]*gz2[k]*0.5*curl_h;
                }

        for(i=MaxX-PML_lx-1;i<MaxX-1;i++)
            for(j=1;j<MaxY;j++)
                for(k=1;k<MaxZ;k++)
                {
                    curl_h=dy_rate*(Hz[i][j][k]-Hz[i][j-1][k])-
                           dz_rate*(Hy[i][j][k]-Hy[i][j][k-1]);
                    Idxh[i-(MaxX-PML_lx-1)][j][k]=
                        Idxh[i-(MaxX-PML_lx-1)][j][k]+curl_h;
                    Dx[i][j][k]=gy3[j]*gz3[k]*Dx[i][j][k]+gy2[j]*gz2[k]*
                        0.5*(curl_h+gx1[i]*Idxh[i-(MaxX-PML_lx-1)][j][k]);
                }

        //Calculate Dy
        for(i=1;i<MaxX;i++)
            for(j=0;j<PML_ly;j++)
                for(k=1;k<MaxZ;k++)
                {
                    curl_h=dz_rate*(Hx[i][j][k]-Hx[i][j][k-1])-
                           dx_rate*(Hz[i][j][k]-Hz[i-1][j][k]);
                    Idyl[i][j][k]=Idyl[i][j][k]+curl_h;
                    Dy[i][j][k]=gx3[i]*gz3[k]*Dy[i][j][k]+gx2[i]*gz2[k]*
                        0.5*(curl_h+gy1[j]*Idyl[i][j][k]);
                }

        for(i=1;i<MaxX;i++)
            for(j=PML_ly;j<MaxY-PML_ly-1;j++)
                for(k=1;k<MaxZ;k++)
                {
                    curl_h=dz_rate*(Hx[i][j][k]-Hx[i][j][k-1])-
                           dx_rate*(Hz[i][j][k]-Hz[i-1][j][k]);
                    Dy[i][j][k]=gx3[i]*gz3[k]*Dy[i][j][k]+
                                gx2[i]*gz2[k]*0.5*curl_h;
                }

        for(i=1;i<MaxX;i++)
            for(j=MaxY-PML_ly-1;j<MaxY-1;j++)
                for(k=1;k<MaxZ;k++)
                {
                    curl_h=dz_rate*(Hx[i][j][k]-Hx[i][j][k-1])-
                           dx_rate*(Hz[i][j][k]-Hz[i-1][j][k]);
                    Idyh[i][j-(MaxY-PML_ly-1)][k]=
                        Idyh[i][j-(MaxY-PML_ly-1)][k]+curl_h;
                    Dy[i][j][k]=gx3[i]*gz3[k]*Dy[i][j][k]+gx2[i]*gz2[k]*
                        0.5*(curl_h+gy1[j]*Idyh[i][j-(MaxY-PML_ly-1)][k]);
                }

        //Calculate Dz
        for(i=1;i<MaxX;i++)
            for(j=1;j<MaxY;j++)
                for(k=0;k<PML_lz;k++)
                {
                    curl_h=dx_rate*(Hy[i][j][k]-Hy[i-1][j][k])-
                           dy_rate*(Hx[i][j][k]-Hx[i][j-1][k]);
                    Idzl[i][j][k]=Idzl[i][j][k]+curl_h;
                    Dz[i][j][k]=gx3[i]*gy3[j]*Dz[i][j][k]+gx2[i]*gy2[j]*
                        0.5*(curl_h+gz1[k]*Idzl[i][j][k]);
                }

        for(i=1;i<MaxX;i++)
            for(j=1;j<MaxY;j++)
                for(k=PML_lz;k<MaxZ-PML_lz-1;k++)
                {
                    curl_h=dx_rate*(Hy[i][j][k]-Hy[i-1][j][k])-
                           dy_rate*(Hx[i][j][k]-Hx[i][j-1][k]);
                    Dz[i][j][k]=gx3[i]*gy3[j]*Dz[i][j][k]+
                                gx2[i]*gy2[j]*0.5*curl_h;
                }

        for(i=1;i<MaxX;i++)
            for(j=1;j<MaxY;j++)
                for(k=MaxZ-PML_lz-1;k<MaxZ-1;k++)
                {
                    curl_h=dx_rate*(Hy[i][j][k]-Hy[i-1][j][k])-
                           dy_rate*(Hx[i][j][k]-Hx[i][j-1][k]);
                    Idzh[i][j][k-(MaxZ-PML_lz-1)]=
                        Idzh[i][j][k-(MaxZ-PML_lz-1)]+curl_h;
                    Dz[i][j][k]=gx3[i]*gy3[j]*Dz[i][j][k]+gx2[i]*gy2[j]*
                        0.5*(curl_h+gz1[k]*Idzh[i][j][k-(MaxZ-PML_lz-1)]);
                }

        return;
}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
//let the out bound of the PML equire zero
void BoundD_2_Zero()
{
    //Normal x (yz plane) , x=0 / x=MaxX-1
    for(j=0;j<MaxY;j++)
        for(k=0;k<MaxZ;k++)
        {
            Dy[0][j][k]=0.0;
            Dz[0][j][k]=0.0;

            Dy[MaxX-1][j][k]=0.0;
            Dz[MaxX-1][j][k]=0.0;
        }

    //Normal y (zx plane) , y=0 / y=MaxY-1
    for(i=0;i<MaxX;i++)
        for(k=0;k<MaxZ;k++)
        {
            Dx[i][0][k]=0.0;
            Dz[i][0][k]=0.0;

            Dx[i][MaxY-1][k]=0.0;
            Dz[i][MaxY-1][k]=0.0;
        }

    //Normal z (xy plane) , z=0 / z=MaxZ-1
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
        {
            Dx[i][j][0]=0.0;
            Dy[i][j][0]=0.0;

            Dx[i][j][MaxZ-1]=0.0;
            Dy[i][j][MaxZ-1]=0.0;
        }

    return;
}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
void AddSourceD()
{
    i=int((SourceX1+SourceX2)/2);
    j=int((SourceY1+SourceY2)/2);
    k=int((SourceZ1+SourceZ2)/2);
    Epsr_Source=Medium_EpsilonR[Model[i][j][k]];

    //Source
        switch(SourceType)
        {
            case 1 :
                SourcePluse=UWB1K*UWB1E0/
                      (exp((UWB1t0-T*dt)/UWB1tr)+exp((T*dt-UWB1t0)/UWB1tp));
                break;
            case 2 :
                SourcePluse=UWB2E0*(exp(-UWB2alpha*T*dt)-exp(-UWB2beita*T*dt));
                break;
            case 3 :
                SourcePluse=exp(-1.*(pow((T*dt-Gausst0)/GaussSpread,2.0)));
                break;
            case 4 :
                SourcePluse=sin(2*Pi*SinFreq*dt*T);
                break;
            case 5:
                SourcePluse=sin(2*Pi*SinFreq*dt*T)*
                            exp(-1.*(pow((T*dt-Gausst0)/GaussSpread,2.0)));
                break;
            case 6:
                SourcePluse=exp(-1.*(Pi*Pi)*(SinFreq*SinFreq)*
                            (T*dt-1./SinFreq)*(T*dt-1./SinFreq));
                SourcePluse=-2.*(Pi*Pi)*(SinFreq*SinFreq)*
                            (T*dt-1./SinFreq)*SourcePluse;
                SourcePluse=SourcePluse*1.E-9;
                break;
            case 7:
                if((T*dt>=RecSint1)&(T*dt<=RecSint2))
                    SourcePluse=sin(2.*Pi*RecSinf0*dt*T);
                else
                    SourcePluse=0.0;
                break;
            default :
                printf("SourceType Error\n");
                return;
        }

        fprintf(fSource,"%16.8E  %16.8E\n",dt*T*1e12,SourcePluse);

        switch(SourceNormal)
        {
            case 3:
                if(SourceVector==1)
                {
                    SourcePluse=SourcePluse/(SourceX2-SourceX1+1)/dx*Epsr_Source;
                    for(i=SourceX1;i<=SourceX2;i++)
                        for(j=SourceY1;j<=SourceY2+1;j++)
                            Dx[i][j][SourceZ1]=Dx[i][j][SourceZ1]+SourcePluse;
                }
                else if(SourceVector==2)
                {
                    SourcePluse=SourcePluse/(SourceY2-SourceY1+1)/dy*Epsr_Source;
                    for(i=SourceX1;i<=SourceX2+1;i++)
                        for(j=SourceY1;j<=SourceY2;j++)
                            Dy[i][j][SourceZ1]=Dy[i][j][SourceZ1]+SourcePluse;
                }
                break;

            case 1:
                if(SourceVector==2)
                {
                    SourcePluse=SourcePluse/(SourceY2-SourceY1+1)/dy*Epsr_Source;
                    for(j=SourceY1;j<=SourceY2;j++)
                        for(k=SourceZ1;k<=SourceZ2+1;k++)
                            Dy[SourceX1][j][k]=Dy[SourceX1][j][k]+SourcePluse;
                }
                else if(SourceVector==3)
                {
                    SourcePluse=SourcePluse/(SourceZ2-SourceZ1+1)/dz*Epsr_Source;
                    for(j=SourceY1;j<=SourceY2+1;j++)
                        for(k=SourceZ1;k<=SourceZ2;k++)
                            Dz[SourceX1][j][k]=Dz[SourceX1][j][k]+SourcePluse;
                }
                break;

            case 2:
                if(SourceVector==1)
                {
                    SourcePluse=SourcePluse/(SourceX2-SourceX1+1)/dx*Epsr_Source;
                    for(i=SourceX1;i<=SourceX2;i++)
                        for(k=SourceZ1;k<=SourceZ2+1;k++)
                            Dx[i][SourceY1][k]=Dx[i][SourceY1][k]+SourcePluse;
                }
                else if(SourceVector==3)
                {
                    SourcePluse=SourcePluse/(SourceZ2-SourceZ1+1)/dz*Epsr_Source;
                    for(i=SourceX1;i<=SourceX2+1;i++)
                        for(k=SourceZ1;k<=SourceZ2;k++)
                            Dz[i][SourceY2][k]=Dz[i][SourceY2][k]+SourcePluse;
                }
                break;

            case 0:
                if(SourceVector==1)
                {
                    for(i=SourceX1;i<=SourceX2;i++)
                        for(k=SourceZ1;k<=SourceZ2;k++)
                            Dx[i][SourceY1][k]=Dx[i][SourceY1][k]+SourcePluse;
                }
                else if(SourceVector==2)
                {
                    for(j=SourceY1;j<=SourceY2;j++)
                        for(k=SourceZ1;k<=SourceZ2;k++)
                            Dy[SourceX1][j][k]=Dy[SourceX1][j][k]+SourcePluse;
                }
                else if(SourceVector==3)
                {
                    for(i=SourceX1;i<=SourceX2;i++)
                        for(k=SourceZ1;k<=SourceZ2;k++)
                            Dz[i][SourceY2][k]=Dz[i][SourceY2][k]+SourcePluse;
                }
                break;

        }

        return;

}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
void AddSoftSource()
{
    i=int((SourceX1+SourceX2)/2);
    j=int((SourceY1+SourceY2)/2);
    k=int((SourceZ1+SourceZ2)/2);
    Epsr_Source=Medium_EpsilonR[Model[i][j][k]];

    //Source
        switch(SourceType)
        {
            case 1 :
                SourcePluse=UWB1K*UWB1E0/
                      (exp((UWB1t0-T*dt)/UWB1tr)+exp((T*dt-UWB1t0)/UWB1tp));
                break;
            case 2 :
                SourcePluse=UWB2E0*(exp(-UWB2alpha*T*dt)-exp(-UWB2beita*T*dt));
                break;
            case 3 :
                SourcePluse=exp(-1.*(pow((T*dt-Gausst0)/GaussSpread,2.0)));
                break;
            case 4 :
                SourcePluse=sin(2*Pi*SinFreq*dt*T);
                break;
            case 5:
                SourcePluse=sin(2*Pi*SinFreq*dt*T)*
                            exp(-1.*(pow((T*dt-Gausst0)/GaussSpread,2.0)));
                break;
            case 6:
                SourcePluse=exp(-1.*(Pi*Pi)*(SinFreq*SinFreq)*
                            (T*dt-1./SinFreq)*(T*dt-1./SinFreq));
                SourcePluse=-2.*(Pi*Pi)*(SinFreq*SinFreq)*
                            (T*dt-1./SinFreq)*SourcePluse;
                SourcePluse=SourcePluse*1.E-9;
                break;
            case 7:
                if((T*dt>=RecSint1)&(T*dt<=RecSint2))
                    SourcePluse=sin(2.*Pi*RecSinf0*dt*T);
                else
                    SourcePluse=0.0;
                break;
            default :
                printf("SourceType Error\n");
                return;
        }

        fprintf(fSource,"%16.8E  %16.8E\n",dt*T*1e12,SourcePluse);

        switch(SourceVector)
        {
			case 1:
                for(i=SourceX1;i<=SourceX2;i++)
                    for(j=SourceY1;j<=SourceY2+1;j++)
						for(k=SourceZ1;k<=SourceZ2;k++)
							Ex[i][j][k]=Ex[i][j][k]+SourcePluse;
				break;

			case 2:
                for(i=SourceX1;i<=SourceX2;i++)
                    for(j=SourceY1;j<=SourceY2+1;j++)
						for(k=SourceZ1;k<=SourceZ2;k++)
							Ey[i][j][k]=Ey[i][j][k]+SourcePluse;
				break;

			case 3:
                for(i=SourceX1;i<=SourceX2;i++)
                    for(j=SourceY1;j<=SourceY2+1;j++)
						for(k=SourceZ1;k<=SourceZ2;k++)
							Ez[i][j][k]=Ez[i][j][k]+SourcePluse;
				break;
        }

        return;

}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
//updata E field
void CalculateE()
{
//Calculate E Field
        ExMax=0.0;
        EyMax=0.0;
        EzMax=0.0;

    for(i=0;i<MaxX-1;i++)
        for(j=0;j<MaxY;j++)
            for(k=0;k<MaxZ;k++)
            {
                Ex[i][j][k]=gax[i][j][k]*
                    (Dx[i][j][k]-Ix[i][j][k]-gcx[i][j][k]*Sx[i][j][k]);
                Ix[i][j][k]=Ix[i][j][k]+gbx[i][j][k]*Ex[i][j][k];
                Sx[i][j][k]=gcx[i][j][k]*Sx[i][j][k]+gdx[i][j][k]*Ex[i][j][k];

                if(ExMax<fabs(Ex[i][j][k]))  ExMax=fabs(Ex[i][j][k]);
            }

    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY-1;j++)
            for(k=0;k<MaxZ;k++)
            {
                Ey[i][j][k]=gay[i][j][k]*
                        (Dy[i][j][k]-Iy[i][j][k]-gcy[i][j][k]*Sy[i][j][k]);
                Iy[i][j][k]=Iy[i][j][k]+gby[i][j][k]*Ey[i][j][k];
                Sy[i][j][k]=gcy[i][j][k]*Sy[i][j][k]+gdy[i][j][k]*Ey[i][j][k];

                if(EyMax<fabs(Ey[i][j][k]))  EyMax=fabs(Ey[i][j][k]);
            }

    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            for(k=0;k<MaxZ-1;k++)
            {
                Ez[i][j][k]=gaz[i][j][k]*
                        (Dz[i][j][k]-Iz[i][j][k]-gcz[i][j][k]*Sz[i][j][k]);
                Iz[i][j][k]=Iz[i][j][k]+gbz[i][j][k]*Ez[i][j][k];
                Sz[i][j][k]=gcz[i][j][k]*Sz[i][j][k]+gdz[i][j][k]*Ez[i][j][k];

                if(EzMax<fabs(Ez[i][j][k]))  EzMax=fabs(Ez[i][j][k]);
            }

    return;

}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
//let electric field in PEC equire zero
void ConductorE_2_Zero()
{
    //Let E on Conductive Element to Zero
    for(k=0;k<MaxZ-1;k++)
        for(j=0;j<MaxY-1;j++)
            for(i=0;i<MaxX-1;i++)
            {
                Sigma=Medium_Sigma[Model[i][j][k]];
                Conform=Medium_Conform[Model[i][j][k]];

                if((Conform==0)&&(Sigma<0.0))
                {
                    Ex[i][j][k]=0.0;
                    Ex[i][j+1][k]=0.0;
                    Ex[i][j][k+1]=0.0;
                    Ex[i][j+1][k+1]=0.0;

                    Ey[i][j][k]=0.0;
                    Ey[i+1][j][k]=0.0;
                    Ey[i][j][k+1]=0.0;
                    Ey[i+1][j][k+1]=0.0;

                    Ez[i][j][k]=0.0;
                    Ez[i+1][j][k]=0.0;
                    Ez[i][j+1][k]=0.0;
                    Ez[i+1][j+1][k]=0.0;
                }

                if(Conform!=0)
                {
                    switch(Conform)
                    {
                        case 1:
                            if(ConType==11)
                            {
                                Ex[i][j][k]=0.0;
                                Ey[i+1][j][k]=0.0;
                            }
                            else if(ConType==12)
                            {
                                Ex[i][j+1][k]=0.0;
                                Ey[i][j][k]=0.0;
                            }
                            else if(ConType==21)
                            {
                                Ex[i][j][k]=0.0;
                                Ey[i][j][k]=0.0;
                            }
                            else if(ConType==22)
                            {
                                Ex[i][j+1][k]=0.0;
                                Ey[i+1][j][k]=0.0;
                            }
                        break;

                        case 2:
                            if(ConType==11)
                            {
                                Ey[i][j][k]=0.0;
                                Ez[i][j+1][k]=0.0;
                            }
                            else if(ConType==12)
                            {
                                Ey[i][j][k+1]=0.0;
                                Ez[i][j][k]=0.0;
                            }
                            else if(ConType==21)
                            {
                                Ey[i][j][k]=0.0;
                                Ez[i][j][k]=0.0;
                            }
                            else if(ConType==22)
                            {
                                Ey[i][j][k+1]=0.0;
                                Ez[i][j+1][k]=0.0;
                            }
                        break;

                        case 3:
                            if(ConType==11)
                            {
                                Ez[i][j][k]=0.0;
                                Ex[i][j][k+1]=0.0;
                            }
                            else if(ConType==12)
                            {
                                Ez[i+1][j][k]=0.0;
                                Ex[i][j][k]=0.0;
                            }
                            else if(ConType==21)
                            {
                                Ez[i][j][k]=0.0;
                                Ex[i][j][k]=0.0;
                            }
                            else if(ConType==22)
                            {
                                Ez[i+1][j][k]=0.0;
                                Ex[i][j][k+1]=0.0;
                            }
                        break;
                    }
                }
            }


    //Let E on Conductive Plate normal x/y/z aixl to Zero
    for(i=0;i<NumNormalX;i++)
    {
        for(j=NormalX[i][1];j<=NormalX[i][2];j++)
            for(k=NormalX[i][3];k<=NormalX[i][4]+1;k++)
                Ey[NormalX[i][0]][j][k]=0.0;


        for(j=NormalX[i][1];j<=NormalX[i][2]+1;j++)
            for(k=NormalX[i][3];k<=NormalX[i][4];k++)
                Ez[NormalX[i][0]][j][k]=0.0;
    }

    for(j=0;j<NumNormalY;j++)
    {
        for(i=NormalY[j][1];i<=NormalY[j][2];i++)
            for(k=NormalY[j][3];k<=NormalY[j][4]+1;k++)
                Ex[i][NormalY[j][0]][k]=0.0;

        for(i=NormalY[j][1];i<=NormalY[j][2]+1;i++)
            for(k=NormalY[j][3]+1;k<=NormalY[j][4];k++)
                Ez[i][NormalY[j][0]][k]=0.0;
    }

    for(k=0;k<NumNormalZ;k++)
    {
        for(i=NormalZ[k][1];i<=NormalZ[k][2];i++)
            for(j=NormalZ[k][3];j<=NormalZ[k][4]+1;j++)
                Ex[i][j][NormalZ[k][0]]=0.0;

        for(i=NormalZ[k][1];i<=NormalZ[k][2]+1;i++)
            for(j=NormalZ[k][3];j<=NormalZ[k][4];j++)
                Ey[i][j][NormalZ[k][0]]=0.0;
    }

    return;

}

//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
//let e field in the out bound of PML equire zero
void BoundE_2_Zero()
{
    //Normal x (yz plane) , x=0 / x=MaxX-1
    for(j=0;j<MaxY-1;j++)
        for(k=0;k<MaxZ;k++)
        {
            Ey[0][j][k]=0.0;
            Ey[MaxX-1][j][k]=0.0;
        }
    for(j=0;j<MaxY;j++)
        for(k=0;k<MaxZ-1;k++)
        {
            Ez[0][j][k]=0.0;
            Ez[MaxX-1][j][k]=0.0;
        }

    //Normal y (zx plane) , y=0 / y=MaxY-1
    for(i=0;i<MaxX-1;i++)
        for(k=0;k<MaxZ;k++)
        {
            Ex[i][0][k]=0.0;
            Ex[i][MaxY-1][k]=0.0;
        }
    for(i=0;i<MaxX;i++)
        for(k=0;k<MaxZ-1;k++)
        {
            Ez[i][0][k]=0.0;
            Ez[i][MaxY-1][k]=0.0;
        }

    //Normal z (xy plane) , z=0 / z=MaxZ-1
    for(i=0;i<MaxX-1;i++)
        for(j=0;j<MaxY;j++)
        {
            Ex[i][j][0]=0.0;
            Ex[i][j][MaxZ-1]=0.0;
        }
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY-1;j++)
        {
            Ey[i][j][0]=0.0;
            Ey[i][j][MaxZ-1]=0.0;
        }

    return;
}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
void AddSource()
{
//Source
        switch(SourceType)
        {
            case 1 :
                SourcePluse=UWB1K*UWB1E0/
                      (exp((UWB1t0-T*dt)/UWB1tr)+exp((T*dt-UWB1t0)/UWB1tp));
                break;
            case 2 :
                SourcePluse=UWB2E0*(exp(-UWB2alpha*T*dt)-exp(-UWB2beita*T*dt));
                break;
            case 3 :
                SourcePluse=exp(-1.*(pow((T*dt-Gausst0)/GaussSpread,2.0)));
                break;
            case 4 :
                SourcePluse=sin(2*Pi*SinFreq*dt*T);
                break;
            case 5:
                SourcePluse=sin(2*Pi*SinFreq*dt*T)*
                            exp(-1.*(pow((T*dt-Gausst0)/GaussSpread,2.0)));
                break;
            case 6:
                SourcePluse=exp(-1.*(Pi*Pi)*(SinFreq*SinFreq)*
                            (T*dt-1./SinFreq)*(T*dt-1./SinFreq));
                SourcePluse=-2.*(Pi*Pi)*(SinFreq*SinFreq)*
                            (T*dt-1./SinFreq)*SourcePluse;
                SourcePluse=SourcePluse*1.E-9;
                break;
            case 7:
                if((T*dt>=RecSint1)&(T*dt<=RecSint2))
                    SourcePluse=sin(2.*Pi*RecSinf0*dt*T);
                else
                    SourcePluse=0.0;
                break;
            default :
                printf("SourceType Error\n");
                return;
        }

        fprintf(fSource,"%16.8E  %16.8E\n",dt*T*1e12,SourcePluse);

        switch(SourceNormal)
        {
            case 3:
                if(SourceVector==1)
                {
                    SourcePluse=SourcePluse/(SourceX2-SourceX1+1)/dx;
                    for(i=SourceX1;i<=SourceX2;i++)
                        for(j=SourceY1;j<=SourceY2+1;j++)
                        {
                            if(SourceHS==0)
                                Ex[i][j][SourceZ1]=SourcePluse;
                            else if(SourceHS==1)
                                Ex[i][j][SourceZ1]=Ex[i][j][SourceZ1]+SourcePluse;
                        }
                }
                else if(SourceVector==2)
                {
                    SourcePluse=SourcePluse/(SourceY2-SourceY1+1)/dy;
                    for(i=SourceX1;i<=SourceX2+1;i++)
                        for(j=SourceY1;j<=SourceY2;j++)
                        {
                            if(SourceHS==0)
                                Ey[i][j][SourceZ1]=SourcePluse;
                            else if(SourceHS==1)
                                Ey[i][j][SourceZ1]=Ey[i][j][SourceZ1]+SourcePluse;
                        }
                }
                break;

            case 1:
                if(SourceVector==2)
                {
                    SourcePluse=SourcePluse/(SourceY2-SourceY1+1)/dy;
                    for(j=SourceY1;j<=SourceY2;j++)
                        for(k=SourceZ1;k<=SourceZ2+1;k++)
                        {
                            if(SourceHS==0)
                                Ey[SourceX1][j][k]=SourcePluse;
                            else if(SourceHS==1)
                                Ey[SourceX1][j][k]=Ey[SourceX1][j][k]+SourcePluse;
                        }
                }
                else if(SourceVector==3)
                {
                    SourcePluse=SourcePluse/(SourceZ2-SourceZ1+1)/dz;
                    for(j=SourceY1;j<=SourceY2+1;j++)
                        for(k=SourceZ1;k<=SourceZ2;k++)
                        {
                            if(SourceHS==0)
                                Ez[SourceX1][j][k]=SourcePluse;
                            else if(SourceHS==1)
                                Ez[SourceX1][j][k]=Ez[SourceX1][j][k]+SourcePluse;
                        }
                }
                break;

            case 2:
                if(SourceVector==1)
                {
                    SourcePluse=SourcePluse/(SourceX2-SourceX1+1)/dx;
                    for(i=SourceX1;i<=SourceX2;i++)
                        for(k=SourceZ1;k<=SourceZ2+1;k++)
                        {
                            if(SourceHS==0)
                                Ex[i][SourceY1][k]=SourcePluse;
                            else if(SourceHS==1)
                                Ex[i][SourceY1][k]=Ex[i][SourceY1][k]+SourcePluse;
                        }
                }
                else if(SourceVector==3)
                {
                    SourcePluse=SourcePluse/(SourceZ2-SourceZ1+1)/dz;
                    for(i=SourceX1;i<=SourceX2+1;i++)
                        for(k=SourceZ1;k<=SourceZ2;k++)
                        {
                            if(SourceHS==0)
                                Ez[i][SourceY1][k]=SourcePluse;
                            else if(SourceHS==1)
                                Ez[i][SourceY2][k]=Ez[i][SourceY2][k]+SourcePluse;
                        }
                }
                break;

            case 0:
                if(SourceVector==1)
                {
                    for(i=SourceX1;i<=SourceX2;i++)
                        for(k=SourceZ1;k<=SourceZ2;k++)
                        {
                            if(SourceHS==0)
                                Ex[i][SourceY1][k]=SourcePluse/dx;
                            else if(SourceHS==1)
                                Ex[i][SourceY1][k]=Ex[i][SourceY1][k]+SourcePluse;
                        }
                }
                else if(SourceVector==2)
                {
                    for(j=SourceY1;j<=SourceY2;j++)
                        for(k=SourceZ1;k<=SourceZ2;k++)
                        {
                            if(SourceHS==0)
                                Ey[SourceX1][j][k]=SourcePluse/dy;
                            else if(SourceHS==1)
                                Ey[SourceX1][j][k]=Ey[SourceX1][j][k]+SourcePluse;
                        }
                }
                else if(SourceVector==3)
                {
                    for(i=SourceX1;i<=SourceX2;i++)
                        for(k=SourceZ1;k<=SourceZ2;k++)
                        {
                            if(SourceHS==0)
                                Ez[i][SourceY1][k]=SourcePluse/dz;
                            else if(SourceHS==1)
                                Ez[i][SourceY2][k]=Ez[i][SourceY2][k]+SourcePluse;
                        }
                }
                break;

        }

        return;

}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
void AddSourceE()
{
    i=int((SourceX1+SourceX2)/2);
    j=int((SourceY1+SourceY2)/2);
    k=int((SourceZ1+SourceZ2)/2);
    Epsr_Source=Medium_EpsilonR[Model[i][j][k]];
    Epsr_Source=1.0;


    switch(SourceNormal)
    {
        case 1:
            if(SourceVector==2)
                for(i=1;i<MaxX;i++)
                    E_inc[i]=gx3[i]*E_inc[i]-
                        gx2[i]*0.5*dx_rate/Epsr_Source*(H_inc[i]-H_inc[i-1]);
            else if(SourceVector==3)
                for(i=1;i<MaxX;i++)
                    E_inc[i]=gx3[i]*E_inc[i]+
                        gx2[i]*0.5*dx_rate/Epsr_Source*(H_inc[i]-H_inc[i-1]);
            break;

        case 2:
            if(SourceVector==3)
                for(j=1;j<MaxY;j++)
                    E_inc[j]=gy3[j]*E_inc[j]-
                        gy2[j]*0.5*dy_rate/Epsr_Source*(H_inc[j]-H_inc[j-1]);
            else if(SourceVector==1)
                for(j=1;j<MaxY;j++)
                    E_inc[j]=gy3[j]*E_inc[j]+
                        gy2[j]*0.5*dy_rate/Epsr_Source*(H_inc[j]-H_inc[j-1]);
            break;

        case 3:
            if(SourceVector==1)
                for(k=1;k<MaxZ;k++)
                    E_inc[k]=gz3[k]*E_inc[k]-
                        gz2[k]*0.5*dz_rate/Epsr_Source*(H_inc[k]-H_inc[k-1]);
            else if(SourceVector==2)
                for(k=1;k<MaxZ;k++)
                    E_inc[k]=gz3[k]*E_inc[k]+
                        gz2[k]*0.5*dz_rate/Epsr_Source*(H_inc[k]-H_inc[k-1]);
            break;
    }

    switch(SourceType)
    {
        case 1 :
            SourcePluse=UWB1K*UWB1E0/
                  (exp((UWB1t0-T*dt)/UWB1tr)+exp((T*dt-UWB1t0)/UWB1tp));
            break;

        case 2 :
            SourcePluse=UWB2E0*(exp(-UWB2alpha*T*dt)-exp(-UWB2beita*T*dt));
            break;

        case 3 :
            SourcePluse=exp(-1.*(pow((T*dt-Gausst0)/GaussSpread,2.0)));
            break;

        case 4 :
            SourcePluse=sin(2*Pi*SinFreq*dt*T);
            break;

        case 5:
            SourcePluse=sin(2*Pi*SinFreq*dt*T)*
                        exp(-1.*(pow((T*dt-Gausst0)/GaussSpread,2.0)));
            break;

        case 6:
            SourcePluse=exp(-1.*(Pi*Pi)*(SinFreq*SinFreq)*
                        (T*dt-1./SinFreq)*(T*dt-1./SinFreq));
            SourcePluse=-2.*(Pi*Pi)*(SinFreq*SinFreq)*
                        (T*dt-1./SinFreq)*SourcePluse;
            SourcePluse=SourcePluse*1.E-9;
            break;

        case 7:
            if((T*dt>=RecSint1)&(T*dt<=RecSint2))
                SourcePluse=sin(2.*Pi*RecSinf0*dt*T);
            else
                SourcePluse=0.0;
            break;

        default :
            printf("SourceType Error\n");
            return;
    }

    fprintf(fSource,"%16.8E  %16.8E\n",dt*T*1e12,SourcePluse);

    switch(SourceNormal)
    {
        case 3:
            if(SourceVector==1)
            {
                SourcePluse=SourcePluse/(SourceX2-SourceX1+1)/dx;
                E_inc[SourceZ1-2]=SourcePluse;
            }
            else if(SourceVector==2)
            {
                SourcePluse=SourcePluse/(SourceY2-SourceY1+1)/dy;
                E_inc[SourceZ1-2]=SourcePluse;
            }
            break;

        case 1:
            if(SourceVector==2)
            {
                SourcePluse=SourcePluse/(SourceY2-SourceY1+1)/dy;
                E_inc[SourceX1-2]=SourcePluse;
            }
            else if(SourceVector==3)
            {
                SourcePluse=SourcePluse/(SourceZ2-SourceZ1+1)/dz;
                E_inc[SourceX1-2]=SourcePluse;
            }
            break;

        case 2:
            if(SourceVector==1)
            {
                SourcePluse=SourcePluse/(SourceX2-SourceX1+1)/dx;
                E_inc[SourceY1-2]=SourcePluse;
            }
            else if(SourceVector==3)
            {
                SourcePluse=SourcePluse/(SourceZ2-SourceZ1+1)/dz;
                E_inc[SourceY1-2]=SourcePluse;
            }
            break;
    }

    switch(SourceNormal)
    {
        case 1:
            if(SourceVector==2)
                for(j=SourceY1;j<=SourceY2;j++)
                    for(k=SourceZ1;k<=SourceZ2+1;k++)
                        Dy[SourceX1][j][k]=Dy[SourceX1][j][k]+
                            0.5*dx_rate/Epsr_Source*H_inc[SourceX1-1];
            else if(SourceVector==3)
                for(j=SourceY1;j<=SourceY2+1;j++)
                    for(k=SourceZ1;k<=SourceZ2;k++)
                        Dz[SourceX1][j][k]=Dz[SourceX1][j][k]-
                            0.5*dx_rate/Epsr_Source*H_inc[SourceX1-1];
            break;

        case 2:
            if(SourceVector==3)
                for(i=SourceX1;i<=SourceX2+1;i++)
                    for(k=SourceZ1;k<=SourceZ2;k++)
                        Dz[i][SourceY1][k]=Dz[i][SourceY1][k]+
                            0.5*dy_rate/Epsr_Source*H_inc[SourceY1-1];
            else if(SourceVector==1)
                for(i=SourceX1;i<=SourceX2;i++)
                    for(k=SourceZ1;k<=SourceZ2+1;k++)
                        Dx[i][SourceY1][k]=Dx[i][SourceY1][k]-
                            0.5*dy_rate/Epsr_Source*H_inc[SourceY1-1];
            break;

        case 3:
            if(SourceVector==1)
                for(i=SourceX1;i<=SourceX2;i++)
                    for(j=SourceY1;j<=SourceY2+1;j++)
                        Dx[i][j][SourceZ1]=Dx[i][j][SourceZ1]+
                            0.5*dz_rate/Epsr_Source*H_inc[SourceZ1-1];
            else if(SourceVector==2)
                for(i=SourceX1;i<=SourceX2+1;i++)
                    for(j=SourceY1;j<=SourceY2;j++)
                        Dy[i][j][SourceZ1]=Dy[i][j][SourceZ1]-
                            0.5*dz_rate/Epsr_Source*H_inc[SourceZ1-1];
            break;
    }

    return;
}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
//save the electric voltage and electric current
//at the appointed ports
void SavePortData()
{
    if(NumPort>0)
    {
    fprintf(fVPort,"%16.8E ",dt*T*1e12);
    fprintf(fIPort1,"%16.8E ",dt*T*1e12);
    fprintf(fIPort2,"%16.8E ",dt*T*1e12);

    for(iPort=0;iPort<NumPort;iPort++)
    {
        Vtemp=0.;
        Itemp1=0.;
        Itemp2=0.;

        if((PortVX1[iPort]!=PortVX2[iPort])&&
           (PortVY1[iPort]==PortVY2[iPort])&&
           (PortVZ1[iPort]==PortVZ2[iPort]))
            PortVector=1;
        else if((PortVX1[iPort]==PortVX2[iPort])&&
                (PortVY1[iPort]!=PortVY2[iPort])&&
                (PortVZ1[iPort]==PortVZ2[iPort]))
            PortVector=2;
        else if((PortVX1[iPort]==PortVX2[iPort])&&
                (PortVY1[iPort]==PortVY2[iPort])&&
                (PortVZ1[iPort]!=PortVZ2[iPort]))
            PortVector=3;
        else
            PortVector=SourceVector;

        if((PortIX1[iPort]!=PortIX2[iPort])&&
           (PortIY1[iPort]!=PortIY2[iPort])&&
           (PortIZ1[iPort]==PortIZ2[iPort]))
            PortNormal=3;
        else if((PortIX1[iPort]==PortIX2[iPort])&&
                (PortIY1[iPort]!=PortIY2[iPort])&&
                (PortIZ1[iPort]!=PortIZ2[iPort]))
            PortNormal=1;
        else if((PortIX1[iPort]!=PortIX2[iPort])&&
                (PortIY1[iPort]==PortIY2[iPort])&&
                (PortIZ1[iPort]!=PortIZ2[iPort]))
            PortNormal=2;
        else
            PortNormal=SourceNormal;


        switch(PortVector)
        {
            case 1:
                for(i=PortVX1[iPort];i<=PortVX2[iPort];i++)
                {
                    ixtemp=i;
                    jytemp=int((PortVY1[iPort]+PortVY2[iPort])/2);
                    kztemp=int((PortVZ1[iPort]+PortVZ2[iPort])/2);
                    Vtemp=Vtemp+Ex[ixtemp][jytemp][kztemp]*dx;
                }
                fprintf(fVPort,"  %16.8E",Vtemp);
                break;

            case 2:
                for(j=PortVY1[iPort];j<=PortVY2[iPort];j++)
                {
                    ixtemp=int((PortVX1[iPort]+PortVX2[iPort])/2);
                    jytemp=j;
                    kztemp=int((PortVZ1[iPort]+PortVZ2[iPort])/2);
                    Vtemp=Vtemp+Ey[ixtemp][jytemp][kztemp]*dy;
                }
                fprintf(fVPort,"  %16.8E",Vtemp);
                break;

            case 3:
                for(k=PortVZ1[iPort];k<=PortVZ2[iPort];k++)
                {
                    ixtemp=int((PortVX1[iPort]+PortVY2[iPort])/2);
                    jytemp=int((PortVY1[iPort]+PortVY2[iPort])/2);
                    kztemp=k;
                    Vtemp=Vtemp+Ez[ixtemp][jytemp][kztemp]*dz;
                }
                fprintf(fVPort,"  %16.8E",Vtemp);
                break;

            default:
                printf("Save VPort error!\n");
                return;
        }

        switch(PortNormal)
        {
            case 3:
                if(PortVector==1)
                {
                    for(j=PortIY1[iPort];j<=PortIY2[iPort]+1;j++)
                    {
                        ixtemp=PortIX1[iPort];
                        kztemp=PortIZ1[iPort];
                        Itemp1=Itemp1-Hy[ixtemp][j][kztemp]*dy;
                        Itemp2=Itemp2-Hy[ixtemp][j][kztemp-1]*dy;
                    }

                    ixtemp=PortIX2[iPort];
                    jytemp=PortIY1[iPort]-1;
                    kztemp=PortIZ1[iPort];
                    Itemp1=Itemp1+Hx[ixtemp][jytemp][kztemp]*dx;
                    Itemp2=Itemp2+Hx[ixtemp][jytemp][kztemp-1]*dx;

                    for(j=PortIY1[iPort];j<=PortIY2[iPort]+1;j++)
                    {
                        ixtemp=PortIX2[iPort];
                        kztemp=PortIZ1[iPort];
                        Itemp1=Itemp1+Hy[ixtemp][j][kztemp]*dy;
                        Itemp2=Itemp2+Hy[ixtemp][j][kztemp-1]*dy;
                    }

                    ixtemp=PortIX2[iPort];
                    jytemp=PortIY2[iPort]+1;
                    kztemp=PortIZ1[iPort];
                    Itemp1=Itemp1-Hx[ixtemp][jytemp][kztemp]*dx;
                    Itemp2=Itemp2-Hx[ixtemp][jytemp][kztemp-1]*dx;

                    Itemp1=-1.*Itemp1;
                    Itemp2=-1.*Itemp2;

                    fprintf(fIPort1,"  %16.8E",Itemp1);
                    fprintf(fIPort2,"  %16.8E",Itemp2);
                }

                else if(PortVector==2)
                {
                    for(i=PortIX1[iPort];i<=PortIX2[iPort]+1;i++)
                    {
                        jytemp=PortIY1[iPort];
                        kztemp=PortIZ1[iPort];
                        Itemp1=Itemp1+Hx[i][jytemp][kztemp]*dx;
                        Itemp2=Itemp2+Hx[i][jytemp][kztemp-1]*dx;
                    }

                    ixtemp=PortIX2[iPort]+1;
                    jytemp=PortIY2[iPort];
                    kztemp=PortIZ1[iPort];
                    Itemp1=Itemp1+Hy[ixtemp][jytemp][kztemp]*dy;
                    Itemp2=Itemp2+Hy[ixtemp][jytemp][kztemp-1]*dy;

                    for(i=PortIX1[iPort];i<=PortIX2[iPort]+1;i++)
                    {
                        jytemp=PortIY2[iPort];
                        kztemp=PortIZ1[iPort];
                        Itemp1=Itemp1-Hx[i][jytemp][kztemp]*dx;
                        Itemp2=Itemp2-Hx[i][jytemp][kztemp-1]*dx;
                    }

                    ixtemp=PortIX1[iPort]-1;
                    jytemp=PortIY2[iPort];
                    kztemp=PortIZ1[iPort];
                    Itemp1=Itemp1-Hy[ixtemp][jytemp][kztemp]*dy;
                    Itemp2=Itemp2-Hy[ixtemp][jytemp][kztemp-1]*dy;

                    Itemp1=-1.*Itemp1;
                    Itemp2=-1.*Itemp2;

                    fprintf(fIPort1,"  %16.8E",Itemp1);
                    fprintf(fIPort2,"  %16.8E",Itemp2);
                }
                break;

            case 1:
                if(PortVector==2)
                {
                    for(k=PortIZ1[iPort];k<=PortIZ2[iPort]+1;k++)
                    {
                        ixtemp=PortIX1[iPort];
                        jytemp=PortIY1[iPort];
                        Itemp1=Itemp1-Hz[ixtemp][jytemp][k]*dz;
                        Itemp2=Itemp2-Hz[ixtemp-1][jytemp][k]*dz;
                    }

                    ixtemp=PortIX1[iPort];
                    jytemp=PortIY2[iPort];
                    kztemp=PortIZ1[iPort]-1;
                    Itemp1=Itemp1+Hy[ixtemp][jytemp][kztemp]*dy;
                    Itemp2=Itemp2+Hy[ixtemp-1][jytemp][kztemp]*dy;

                    for(k=PortIZ1[iPort];k<=PortIZ2[iPort];k++)
                    {
                        ixtemp=PortIX1[iPort];
                        jytemp=PortIY2[iPort];
                        Itemp1=Itemp1+Hz[ixtemp][jytemp][k]*dz;
                        Itemp2=Itemp2+Hz[ixtemp-1][jytemp][k]*dz;
                    }

                    ixtemp=PortIX1[iPort];
                    jytemp=PortIY2[iPort];
                    kztemp=PortIZ2[iPort]+1;
                    Itemp1=Itemp1-Hy[ixtemp][jytemp][kztemp]*dy;
                    Itemp2=Itemp2-Hy[ixtemp-1][jytemp][kztemp]*dy;

                    Itemp1=-1.*Itemp1;
                    Itemp2=-1.*Itemp2;

                    fprintf(fIPort1,"  %16.8E",Itemp1);
                    fprintf(fIPort2,"  %16.8E",Itemp2);
                }

                else if(PortVector==3)
                {
                    for(j=PortIY1[iPort];j<=PortIY2[iPort]+1;j++)
                    {
                        ixtemp=PortIX1[iPort];
                        kztemp=PortIZ1[iPort];
                        Itemp1=Itemp1+Hy[ixtemp][j][kztemp]*dy;
                        Itemp2=Itemp2+Hy[ixtemp-1][j][kztemp]*dy;
                    }

                    ixtemp=PortIX1[iPort];
                    jytemp=PortIY1[iPort]-1;
                    kztemp=PortIZ2[iPort];
                    Itemp1=Itemp1-Hz[ixtemp][jytemp][kztemp]*dz;
                    Itemp2=Itemp2-Hz[ixtemp-1][jytemp][kztemp]*dz;

                    for(j=PortIY1[iPort];j<=PortIY2[iPort]+1;j++)
                    {
                        ixtemp=PortIX1[iPort];
                        kztemp=PortIZ1[iPort];
                        Itemp1=Itemp1-Hy[ixtemp][j][kztemp]*dy;
                        Itemp2=Itemp2-Hy[ixtemp-1][j][kztemp]*dy;
                    }

                    ixtemp=PortIX1[iPort];
                    jytemp=PortIY2[iPort]+1;
                    kztemp=PortIZ2[iPort];
                    Itemp1=Itemp1+Hz[ixtemp][jytemp][kztemp]*dz;
                    Itemp2=Itemp2+Hz[ixtemp-1][jytemp][kztemp]*dz;

                    Itemp1=-1.*Itemp1;
                    Itemp2=-1.*Itemp2;

                    fprintf(fIPort1,"  %16.8E",Itemp1);
                    fprintf(fIPort2,"  %16.8E",Itemp2);
                }
                break;

            case 2:
                if(PortVector==1)
                {
                    for(k=PortIZ1[iPort];k<=PortIZ2[iPort]+1;k++)
                    {
                        ixtemp=PortIX1[iPort];
                        jytemp=PortIY1[iPort];
                        Itemp1=Itemp1+Hz[ixtemp][jytemp][k]*dz;
                        Itemp2=Itemp2+Hz[ixtemp][jytemp-1][k]*dz;
                    }

                    ixtemp=PortIX2[iPort];
                    jytemp=PortIY1[iPort];
                    kztemp=PortIZ1[iPort]-1;
                    Itemp1=Itemp1-Hx[ixtemp][jytemp][kztemp]*dx;
                    Itemp2=Itemp2-Hx[ixtemp][jytemp-1][kztemp]*dx;

                    for(k=PortIZ1[iPort];k<=PortIZ2[iPort]+1;k++)
                    {
                        ixtemp=PortIX2[iPort];
                        jytemp=PortIY1[iPort];
                        Itemp1=Itemp1-Hz[ixtemp][jytemp][k]*dz;
                        Itemp2=Itemp2-Hz[ixtemp][jytemp-1][k]*dz;
                    }

                    ixtemp=PortIX2[iPort];
                    jytemp=PortIY1[iPort];
                    kztemp=PortIZ2[iPort]+1;
                    Itemp1=Itemp1+Hx[ixtemp][jytemp][kztemp]*dx;
                    Itemp2=Itemp2+Hx[ixtemp][jytemp-1][kztemp]*dx;

                    Itemp1=-1.*Itemp1;
                    Itemp2=-1.*Itemp2;

                    fprintf(fIPort1,"  %16.8E",Itemp1);
                    fprintf(fIPort2,"  %16.8E",Itemp2);
                }

                else if(PortVector==3)
                {
                    for(i=PortIX1[iPort];i<=PortIX2[iPort]+1;i++)
                    {
                        jytemp=PortIY1[iPort];
                        kztemp=PortIZ1[iPort];
                        Itemp1=Itemp1-Hx[i][jytemp][kztemp]*dx;
                        Itemp2=Itemp2-Hx[i][jytemp-1][kztemp]*dx;
                    }

                    ixtemp=PortIX1[iPort]-1;
                    jytemp=PortIY1[iPort];
                    kztemp=PortIZ2[iPort];
                    Itemp1=Itemp1+Hz[ixtemp][jytemp][kztemp]*dz;
                    Itemp2=Itemp2+Hz[ixtemp][jytemp-1][kztemp]*dz;

                    for(i=PortIX1[iPort];i<=PortIX2[iPort]+1;i++)
                    {
                        jytemp=PortIY1[iPort];
                        kztemp=PortIZ2[iPort];
                        Itemp1=Itemp1+Hx[i][jytemp][kztemp]*dx;
                        Itemp2=Itemp2+Hx[i][jytemp-1][kztemp]*dx;
                    }

                    ixtemp=PortIX2[iPort]+1;
                    jytemp=PortIY1[iPort];
                    kztemp=PortIZ2[iPort];
                    Itemp1=Itemp1-Hz[ixtemp][jytemp][kztemp]*dz;
                    Itemp2=Itemp2-Hz[ixtemp][jytemp-1][kztemp]*dz;

                    Itemp1=-1.*Itemp1;
                    Itemp2=-1.*Itemp2;

                    fprintf(fIPort1,"  %16.8E",Itemp1);
                    fprintf(fIPort2,"  %16.8E",Itemp2);
                }
                break;

            default:
                printf("Save IPort error!\n");
                return;
        }
    }

    fprintf(fVPort,"\n");
    fprintf(fIPort1,"\n");
    fprintf(fIPort2,"\n");
    }

    return;
}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
//save the electric field in appointed plane
//at appointed time step
void SaveFieldonPlane()
{
    for(iFP=0;iFP<NumFieldPlane;iFP++)
      if(T==FieldPlaneTimeStep[iFP])
          switch(FieldPlaneNormal[iFP])
          {
             case 1:
                 for(j=0;j<MaxY;j++)
                   for(k=0;k<MaxZ;k++)
                     switch(SourceVector)
                     {
                        case 1:
                            fprintf(fFieldPlane[iFP], "%d  %d  %16.8E\n",
                                j,k,Ex[FieldPlaneVector[iFP]][j][k]);
                            break;

                        case 2:
                            fprintf(fFieldPlane[iFP], "%d  %d  %16.8E\n",
                                j,k,Ey[FieldPlaneVector[iFP]][j][k]);
                            break;

                        case 3:
                            fprintf(fFieldPlane[iFP], "%d  %d  %16.8E\n",
                                j,k,Ez[FieldPlaneVector[iFP]][j][k]);
                            break;
                     }
                 break;

             case 2:
                 for(k=0;k<MaxZ;k++)
                   for(i=0;i<MaxX;i++)
                     switch(SourceVector)
                     {
                        case 1:
                            fprintf(fFieldPlane[iFP], "%d  %d  %16.8E\n",
                                k,i,Ex[i][FieldPlaneVector[iFP]][k]);
                            break;

                        case 2:
                            fprintf(fFieldPlane[iFP], "%d  %d  %16.8E\n",
                                k,i,Ey[i][FieldPlaneVector[iFP]][k]);
                            break;

                        case 3:
                            fprintf(fFieldPlane[iFP], "%d  %d  %16.8E\n",
                                k,i,Ez[i][FieldPlaneVector[iFP]][k]);
                            break;
                     }
                 break;

             case 3:
                 for(i=0;i<MaxX;i++)
                   for(j=0;j<MaxY;j++)
                     switch(SourceVector)
                     {
                        case 1:
                            fprintf(fFieldPlane[iFP], "%d  %d  %16.8E\n",
                                i,j,Ex[i][j][FieldPlaneVector[iFP]]);
                            break;

                        case 2:
                            fprintf(fFieldPlane[iFP], "%d  %d  %16.8E\n",
                                i,j,Ey[i][j][FieldPlaneVector[iFP]]);
                            break;

                        case 3:
                            fprintf(fFieldPlane[iFP], "%d  %d  %16.8E\n",
                                i,j,Ez[i][j][FieldPlaneVector[iFP]]);
                            break;
                     }
                 break;

          }
      return;
}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
//updata H field
void CalculateH()
{
        //Calculate Hx
        for(i=0;i<PML_lx;i++)
            for(j=0;j<MaxY-1;j++)
                for(k=0;k<MaxZ-1;k++)
                {
                    curl_e=dz_rate*(Ey[i][j][k+1]-Ey[i][j][k])-
                           dy_rate*(Ez[i][j+1][k]-Ez[i][j][k]);
                    Ihxl[i][j][k]=Ihxl[i][j][k]+curl_e;
                    Hx[i][j][k]=fy3[j]*fz3[k]*Hx[i][j][k]+fy2[j]*fz2[k]*
                        0.5*(curl_e+fx1[i]*Ihxl[i][j][k]);
                }

        for(i=PML_lx;i<MaxX-PML_lx;i++)
            for(j=0;j<MaxY-1;j++)
                for(k=0;k<MaxZ-1;k++)
                {
                    curl_e=dz_rate*(Ey[i][j][k+1]-Ey[i][j][k])-
                           dy_rate*(Ez[i][j+1][k]-Ez[i][j][k]);
                    Hx[i][j][k]=fy3[j]*fz3[k]*Hx[i][j][k]+
                                fy2[j]*fz2[k]*0.5*curl_e;
                }

        for(i=MaxX-PML_lx;i<MaxX;i++)
            for(j=0;j<MaxY-1;j++)
                for(k=0;k<MaxZ-1;k++)
                {
                    curl_e=dz_rate*(Ey[i][j][k+1]-Ey[i][j][k])-
                           dy_rate*(Ez[i][j+1][k]-Ez[i][j][k]);
                    Ihxh[i-MaxX+PML_lx][j][k]=Ihxh[i-MaxX+PML_lx][j][k]+curl_e;
                    Hx[i][j][k]=fy3[j]*fz3[k]*Hx[i][j][k]+fy2[j]*fz2[k]*
                        0.5*(curl_e+fx1[i]*Ihxh[i-MaxX+PML_lx][j][k]);
                }

        //Calculate Hy
        for(i=0;i<MaxX-1;i++)
            for(j=0;j<PML_ly;j++)
                for(k=0;k<MaxZ-1;k++)
                {
                    curl_e=dx_rate*(Ez[i+1][j][k]-Ez[i][j][k])-
                           dz_rate*(Ex[i][j][k+1]-Ex[i][j][k]);
                    Ihyl[i][j][k]=Ihyl[i][j][k]+curl_e;
                    Hy[i][j][k]=fx3[i]*fz3[k]*Hy[i][j][k]+fx2[i]*fz2[k]*
                        0.5*(curl_e+fy1[j]*Ihyl[i][j][k]);
                }

        for(i=0;i<MaxX-1;i++)
            for(j=PML_ly;j<MaxY-PML_ly;j++)
                for(k=0;k<MaxZ-1;k++)
                {
                    curl_e=dx_rate*(Ez[i+1][j][k]-Ez[i][j][k])-
                           dz_rate*(Ex[i][j][k+1]-Ex[i][j][k]);
                    Hy[i][j][k]=fx3[i]*fz3[k]*Hy[i][j][k]+
                                fx2[i]*fz2[k]*0.5*curl_e;
                }

        for(i=0;i<MaxX-1;i++)
            for(j=MaxY-PML_ly;j<MaxY;j++)
                for(k=0;k<MaxZ-1;k++)
                {
                    curl_e=dx_rate*(Ez[i+1][j][k]-Ez[i][j][k])-
                           dz_rate*(Ex[i][j][k+1]-Ex[i][j][k]);
                    Ihyh[i][j-MaxY+PML_ly][k]=Ihyh[i][j-MaxY+PML_ly][k]+curl_e;
                    Hy[i][j][k]=fx3[i]*fz3[k]*Hy[i][j][k]+fx2[i]*fz2[k]*
                        0.5*(curl_e+fy1[j]*Ihyh[i][j-MaxY+PML_ly][k]);
                }

        //Calculate Hz
        for(i=0;i<MaxX-1;i++)
            for(j=0;j<MaxY-1;j++)
                for(k=0;k<PML_lz;k++)
                {
                    curl_e=dy_rate*(Ex[i][j+1][k]-Ex[i][j][k])-
                           dx_rate*(Ey[i+1][j][k]-Ey[i][j][k]);
                    Ihzl[i][j][k]=Ihzl[i][j][k]+curl_e;
                    Hz[i][j][k]=fx3[i]*fy3[j]*Hz[i][j][k]+fx2[i]*fy2[j]*
                        0.5*(curl_e+fz1[k]*Ihzl[i][j][k]);
                }

        for(i=0;i<MaxX-1;i++)
            for(j=0;j<MaxY-1;j++)
                for(k=PML_lz;k<MaxZ-PML_lz;k++)
                {
                    curl_e=dy_rate*(Ex[i][j+1][k]-Ex[i][j][k])-
                           dx_rate*(Ey[i+1][j][k]-Ey[i][j][k]);
                    Hz[i][j][k]=fx3[i]*fy3[j]*Hz[i][j][k]+
                                fx2[i]*fy2[j]*0.5*curl_e;
                }

        for(i=0;i<MaxX-1;i++)
            for(j=0;j<MaxY-1;j++)
                for(k=MaxZ-PML_lz;k<MaxZ;k++)
                {
                    curl_e=dy_rate*(Ex[i][j+1][k]-Ex[i][j][k])-
                           dx_rate*(Ey[i+1][j][k]-Ey[i][j][k]);
                    Ihzh[i][j][k-MaxZ+PML_lz]=Ihzh[i][j][k-MaxZ+PML_lz]+curl_e;
                    Hz[i][j][k]=fx3[i]*fy3[j]*Hz[i][j][k]+fx2[i]*fy2[j]*
                        0.5*(curl_e+fz1[k]*Ihzh[i][j][k-MaxZ+PML_lz]);
                }

        return;

}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
void AddSourceH()
{
    i=int((SourceX1+SourceX2)/2);
    j=int((SourceY1+SourceY2)/2);
    k=int((SourceZ1+SourceZ2)/2);
    Epsr_Source=Medium_EpsilonR[Model[i][j][k]];

    switch(SourceNormal)
    {
        case 1:
            if(SourceVector==2)
                for(i=1;i<MaxX;i++)
                    H_inc[i]=fx3[i]*H_inc[i]-
                        fx2[i]*0.5*dx_rate*(E_inc[i+1]-E_inc[i])/Epsr_Source;
            else if(SourceVector==3)
                for(i=1;i<MaxX;i++)
                    H_inc[i]=fx3[i]*H_inc[i]+
                        fx2[i]*0.5*dx_rate*(E_inc[i+1]-E_inc[i])/Epsr_Source;
            break;

        case 2:
            if(SourceVector==3)
                for(j=1;j<MaxY;j++)
                    H_inc[j]=fy3[j]*H_inc[j]-
                        fy2[j]*0.5*dy_rate*(E_inc[j+1]-E_inc[j])/Epsr_Source;
            else if(SourceVector==1)
                for(j=1;j<MaxY;j++)
                    H_inc[j]=fy3[j]*H_inc[j]+
                        fy2[j]*0.5*dy_rate*(E_inc[j+1]-E_inc[j])/Epsr_Source;
            break;

        case 3:
            if(SourceVector==1)
                for(k=1;k<MaxZ;k++)
                    H_inc[k]=fz3[k]*H_inc[k]-
                        fz2[k]*0.5*dz_rate*(E_inc[k+1]-E_inc[k])/Epsr_Source;
            else if(SourceVector==2)
                for(k=1;k<MaxZ;k++)
                    H_inc[k]=fz3[k]*H_inc[k]+
                        fz2[k]*0.5*dz_rate*(E_inc[k+1]-E_inc[k])/Epsr_Source;
            break;
    }

    switch(SourceNormal)
    {
        case 1:
            if(SourceVector==2)
                for(j=SourceY1;j<=SourceY2;j++)
                    for(k=SourceZ1;k<=SourceZ2+1;k++)
                        Hz[SourceX1-1][j][k]=Hz[SourceX1-1][j][k]+
                            0.5*dx_rate*E_inc[SourceX1];
            else if(SourceVector==3)
                for(j=SourceY1;j<=SourceY2+1;j++)
                    for(k=SourceZ1;k<=SourceZ2;k++)
                        Hy[SourceX1-1][j][k]=Hy[SourceX1-1][j][k]-
                            0.5*dx_rate*E_inc[SourceX1];
            break;

        case 2:
            if(SourceVector==3)
                for(i=SourceX1;i<=SourceX2+1;i++)
                    for(k=SourceZ1;k<=SourceZ2;k++)
                        Hx[i][SourceY1-1][k]=Hx[i][SourceY1-1][k]+
                            0.5*dy_rate*E_inc[SourceY1];
            else if(SourceVector==1)
                for(i=SourceX1;i<=SourceX2;i++)
                    for(k=SourceZ1;k<=SourceZ2+1;k++)
                        Hz[i][SourceY1-1][k]=Hz[i][SourceY1-1][k]-
                            0.5*dy_rate*E_inc[SourceY1];
            break;

        case 3:
            if(SourceVector==1)
                for(i=SourceX1;i<=SourceX2;i++)
                    for(j=SourceY1;j<=SourceY2+1;j++)
                        Hy[i][j][SourceZ1-1]=Hy[i][j][SourceZ1-1]+
                            0.5*dz_rate*E_inc[SourceZ1];
            else if(SourceVector==2)
                for(i=SourceX1;i<=SourceX2+1;i++)
                    for(j=SourceY1;j<=SourceY2;j++)
                        Hx[i][j][SourceZ1-1]=Hx[i][j][SourceZ1-1]-
                            0.5*dz_rate*E_inc[SourceZ1];
            break;
    }

    return;
}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
void BoundH_2_Zero()
{
    //Normal x (yz plane) , x=0 / x=MaxX-1
    for(j=0;j<MaxY;j++)
        for(k=0;k<MaxZ;k++)
        {
            Hx[0][j][k]=0.0;
            Hx[MaxX-1][j][k]=0.0;
        }

    //Normal y (zx plane) , y=0 / y=MaxY-1
    for(i=0;i<MaxX;i++)
        for(k=0;k<MaxZ;k++)
        {
            Hy[i][0][k]=0.0;
            Hy[i][MaxY-1][k]=0.0;
        }

    //Normal z (xy plane) , z=0 / z=MaxZ-1
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
        {
            Hz[i][j][0]=0.0;
            Hz[i][j][MaxZ-1]=0.0;
        }

    return;
}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
//clear the system resource the code used, such as memory and file handle.
void DelData()
{
    delete WatchX;
    delete WatchY;
    delete WatchZ;

    for(i=0;i<NumNormalX;i++) delete NormalX[i];
    delete NormalX;
    for(i=0;i<NumNormalY;i++) delete NormalY[i];
    delete NormalY;
    for(i=0;i<NumNormalZ;i++) delete NormalZ[i];
    delete NormalZ;

    //delete Model
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Model[i][j];
    for(i=0;i<MaxX;i++) delete Model[i];
    delete Model;

    //delete Dx,Dy,Dz
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Dx[i][j];
    for(i=0;i<MaxX;i++) delete Dx[i];
    delete Dx;
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Dy[i][j];
    for(i=0;i<MaxX;i++) delete Dy[i];
    delete Dy;
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Dz[i][j];
    for(i=0;i<MaxX;i++) delete Dz[i];
    delete Dz;

    //delete Ex,Ey,Ez
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Ex[i][j];
    for(i=0;i<MaxX;i++) delete Ex[i];
    delete Ex;
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Ey[i][j];
    for(i=0;i<MaxX;i++) delete Ey[i];
    delete Ey;
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Ez[i][j];
    for(i=0;i<MaxX;i++) delete Ez[i];
    delete Ez;

    //delete gax,gay,gaz
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete gax[i][j];
    for(i=0;i<MaxX;i++) delete gax[i];
    delete gax;
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete gay[i][j];
    for(i=0;i<MaxX;i++) delete gay[i];
    delete gay;
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete gaz[i][j];
    for(i=0;i<MaxX;i++) delete gaz[i];
    delete gaz;

    //delete gbx,gby,gbz
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete gbx[i][j];
    for(i=0;i<MaxX;i++) delete gbx[i];
    delete gbx;
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete gby[i][j];
    for(i=0;i<MaxX;i++) delete gby[i];
    delete gby;
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete gbz[i][j];
    for(i=0;i<MaxX;i++) delete gbz[i];
    delete gbz;

    //delete Ix,Iy,Iz
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Ix[i][j];
    for(i=0;i<MaxX;i++) delete Ix[i];
    delete Ix;
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Iy[i][j];
    for(i=0;i<MaxX;i++) delete Iy[i];
    delete Iy;
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Iz[i][j];
    for(i=0;i<MaxX;i++) delete Iz[i];
    delete Iz;

    //delete Hx,Hy,Hz
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Hx[i][j];
    for(i=0;i<MaxX;i++) delete Hx[i];
    delete Hx;
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Hy[i][j];
    for(i=0;i<MaxX;i++) delete Hy[i];
    delete Hy;
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Hz[i][j];
    for(i=0;i<MaxX;i++) delete Hz[i];
    delete Hz;

    //delete HB
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete HB[i][j];
    for(i=0;i<MaxX;i++) delete HB[i];
    delete HB;

    //delete Id
    for(i=0;i<PML_lx;i++)
        for(j=0;j<MaxY;j++)
            delete Idxl[i][j];
    for(i=0;i<PML_lx;i++) delete Idxl[i];
    delete Idxl;
    for(i=0;i<PML_lx;i++)
        for(j=0;j<MaxY;j++)
            delete Idxh[i][j];
    for(i=0;i<PML_lx;i++) delete Idxh[i];
    delete Idxh;

    for(i=0;i<MaxX;i++)
        for(j=0;j<PML_ly;j++)
            delete Idyl[i][j];
    for(i=0;i<MaxX;i++) delete Idyl[i];
    delete Idyl;
    for(i=0;i<MaxX;i++)
        for(j=0;j<PML_ly;j++)
            delete Idyh[i][j];
    for(i=0;i<MaxX;i++) delete Idyh[i];
    delete Idyh;

    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Idzl[i][j];
    for(i=0;i<MaxX;i++) delete Idzl[i];
    delete Idzl;
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Idzh[i][j];
    for(i=0;i<MaxX;i++) delete Idzh[i];
    delete Idzh;

    //delete Ih
    for(i=0;i<PML_lx;i++)
        for(j=0;j<MaxY;j++)
            delete Ihxl[i][j];
    for(i=0;i<PML_lx;i++) delete Ihxl[i];
    delete Ihxl;
    for(i=0;i<PML_lx;i++)
        for(j=0;j<MaxY;j++)
            delete Ihxh[i][j];
    for(i=0;i<PML_lx;i++) delete Ihxh[i];
    delete Ihxh;

    for(i=0;i<MaxX;i++)
        for(j=0;j<PML_ly;j++)
            delete Ihyl[i][j];
    for(i=0;i<MaxX;i++) delete Ihyl[i];
    delete Ihyl;
    for(i=0;i<MaxX;i++)
        for(j=0;j<PML_ly;j++)
            delete Ihyh[i][j];
    for(i=0;i<MaxX;i++) delete Ihyh[i];
    delete Ihyh;

    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Ihzl[i][j];
    for(i=0;i<MaxX;i++) delete Ihzl[i];
    delete Ihzl;
    for(i=0;i<MaxX;i++)
        for(j=0;j<MaxY;j++)
            delete Ihzh[i][j];
    for(i=0;i<MaxX;i++) delete Ihzh[i];
    delete Ihzh;

    //delete gx,gy,gz,fx,fy,fz
    delete gx3;
    delete gx2;
    delete gx1;
    delete fx3;
    delete fx2;
    delete fx1;

    delete gy3;
    delete gy2;
    delete gy1;
    delete fy3;
    delete fy2;
    delete fy1;

    delete gz3;
    delete gz2;
    delete gz1;
    delete fz3;
    delete fz2;
    delete fz1;

    //close file
    fclose(fSource);
	if(NumPort>0)
	{
		fclose(fVPort);
		fclose(fIPort1);
		fclose(fIPort2);
	}
    for(i=0;i<NumFieldPlane;i++)
        fclose(fFieldPlane[i]);

    DelDataFlag=1;
    ReadDataFlag=0;

    return;
}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
void PostProcessing()
{
    PostCalculateZ();

    PostCalculateS();

    PostCalculateBeita();

    return;
}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
void PostCalculateZ()
{

    return;
}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
void PostCalculateS()
{
    return;
}
//---------------------------------------------------------------------------


//---------------------------------------------------------------------------
void PostCalculateBeita()
{
    return;
}
//---------------------------------------------------------------------------


